Environmental DNA (eDNA) metabarcoding differentiates between micro-habitats within the rocky intertidal

### Loading Packages###

renv::restore()

packages <- c(
  "tidyverse",
  "sf",
  "patchwork",
  "stringr",
  "eulerr",
  "indicspecies",
  "cluster",
  "vegan",
  "chisq.posthoc.test",
  "taxize",
  "Polychrome",
  "pals",
  "spocc",
  "mregions",
  "rgbif",
  "mapview",
  "ggbeeswarm",
  "ggh4x",
  "chunkhooks",
  "devtools",
  "grid",
  "styler",
  "phyloseq",
  "ggtree",
  "ggtreeExtra",
  "ampvis2",
  "pdftools",
  "gt",
  "gtExtras"
)

invisible(lapply(packages, library, character.only = TRUE))

hook_figure_unit()

### Core Functions###
# Collection of functions used across analyses

# Function that removes NCBI Taxonomy IDs from the end of taxonomic names (as in the Anacapa Pipeline output)
remove_NCBInums <- function(list) {
  list <- sub("_[^_]+$", "", list)
  return(list)
}

# Function that converts an ampvis2 object to a phyloseq object
# Can be used either on an existing ampvis2 object (abund, tax, md = NULL) or separate abund, tax, and metadata data frames (ampvis = NULL)
amp_to_phyloseq <- function(ampvis, abund, tax, md) {
  if (!is.null(ampvis)) {
    physeq <- phyloseq(
      otu_table(as.matrix(ampvis$abund), taxa_are_rows = TRUE),
      tax_table(as.matrix(ampvis$tax)),
      sample_data(ampvis$metadata)
    )
  } else if (is.null(md)) {
    physeq <- phyloseq(
      otu_table(as.matrix(abund), taxa_are_rows = TRUE),
      tax_table(as.matrix(tax))
    )
  } else {
    physeq <- phyloseq(
      otu_table(as.matrix(abund), taxa_are_rows = TRUE),
      tax_table(as.matrix(tax)),
      sample_data(md)
    )
  }
  return(physeq)
}

# Function that saves .png and .pdf versions of a given figure
save_pdf_png <- function(plot, path, w, h, u, bg, dpi) {
  for (x in c("pdf", "png")) {
    ggsave(
      plot = plot,
      filename = paste0(path, x),
      device = x,
      bg = bg,
      width = w,
      height = h,
      units = u,
      dpi = dpi
    )
  }
}

Materials and Methods

2.1 Detailed Protocols

To improve reproducibility (e.g. Dickie et al., 2018; Shea et al., 2023), enable open data science (e.g. Fredston and Lowndes, 2024), and aid in the initiation of new eDNA biomonitoring projects, we have published detailed, step-by-step protocols for many of the methods described, including specific materials used, photographs, and additional methodological notes not possible to include here. See Shea (2023) for sample collection and filtering, Shea (2023) for DNA extractions, Shea (2023) for PCR amplification, and Shea (2023) for shipping samples. Additionally, we have published all data and code for replicating our analyses via Dryad, including FASTQ files, our modified Anacapa Container and scripts for bioinformatics, unprocessed and processed datasets, and an R Markdown file that reproduces all methods & results detailed here (Dryad citation).

2.2 Sampling Site

To better understand spatial and temporal differences in eDNA signals in a complex coastal environment, we sought a rocky intertidal field location that had consistent, large, accessible tide pools that were fully isolated from one another at some low tides but interconnected during other parts of their exposure period. We selected the intertidal at Pillar Point, a headlands promontory to the west of Pillar Point Harbor in San Mateo County, California, USA. Pillar Point is a popular recreational intertidal site that is directly adjacent to the Pillar Point State Marine Conservation Area (no specific use scientific collection permit required).

Within Pillar Point, we sampled at three discrete locations: two individual tide pools with a range of physical connectivity (Tide Pool 1, S1: 37.495342°, -122.498731°; Tide Pool 2, S2: 37.495000°, -122.498947°) and an equidistant location (Nearshore, N: 37.4952833°, -122.4992056°) where there was well-mixed offshore water for the duration of the tidal cycle (Figure 2; not produced in R). Each site was approximately 40 meters from all other sites. Tide Pool 1 and Tide Pool 2 are fully isolated at tidal heights of around 0 m (mean low low water, MLLW) or lower, and substantially connected at around 0.25 m or higher. On the day we sampled, this meant water actively flowed between the locations at the start (11:30 PST) and end (17:00 PST) of the sampling period, but that the sites were disconnected at low tide in the middle of the sampling period.

2.3 Sample Collection & Filtration

We collected 1 L surface samples from each location every 30 minutes for the duration of time the rocky intertidal was exposed on 28 January 2022, using single-use enteral feeding pouches (Covidien, Dublin, Ireland). Sampling commenced at 11:30 PST; the exact times of sample collection, in relation to the tide, are shown in Figure 3. Following the approach used by Gold et al. (2021b), we attached a sterile 0.22 \(\mu\)m pore size Sterivex cartridge (MilliporeSigma, Burlington, MA, USA) to the tubing of each feeding pouch, allowing samples to be immediately gravity filtered in the field. While gravity filtering (1-2 hours per sample), samples were shaded with an awning to prevent any degradation by sunlight (Andruszkiewicz et al., 2017).

### Figure 3###
# Produce the top portion of Figure 3 from tide data

# Read Pillar Point NOAA/NOS/CO-OPs Tide Chart data file
Tides <- read.table("Data/28-Jan-2022-TideChart.txt", header = TRUE, sep = "", dec = ".", skip = 13)

# Convert date format in table using lubridate
Tides$datetime <- lubridate::ymd_hm(paste(Tides$Date, Tides$Time), tz = "US/Pacific")

# Set bounds of sampling period to highlight
xmin <- Tides$datetime[47]
xmax <- Tides$datetime[69]
ymin <- -.5
ymax <- .5

box_color <- "#808080"

# Create plot of full tidal cycle with sampling period highlighted
p1 <- ggplot(data = Tides, aes(x = datetime, y = Pred, group = 1)) +
  annotate("rect", xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = box_color) +
  geom_line() +
  coord_cartesian(xlim = c(Tides$datetime[1], Tides$datetime[96])) +
  theme_classic(base_size = 12) +
  theme(
    panel.background = element_rect(fill = "transparent", colour = NA_character_),
    # necessary to avoid drawing panel outline
    plot.background = element_rect(fill = "transparent", colour = NA_character_),
    # necessary to avoid drawing plot outline
    axis.text.x = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    axis.text.y = element_text(color = "black"),
    panel.border = element_rect(colour = "black", fill = NA, linewidth = 1)
  ) +
  xlab("Time (PST)") +
  ylab("MLLW Tide Prediction (m)") +
  scale_x_datetime(date_breaks = "2 hour", date_labels = "%H:00", expand = c(0, 0), position = "top")

# Create plot of just sampling period
p2 <- ggplot(data = Tides, aes(x = datetime, y = Pred, group = 1)) +
  geom_line() +
  coord_cartesian(xlim = c(xmin, xmax), ylim = c(ymin, ymax)) +
  theme_classic(base_size = 12) +
  theme(
    panel.background = element_rect(fill = box_color),
    plot.background = element_rect(fill = "transparent", colour = NA_character_),
    axis.text.x = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    axis.text.y = element_text(color = "black"),
    panel.border = element_rect(colour = "black", fill = NA, linewidth = 1)
  ) +
  theme(axis.title.x = element_blank()) +
  ylab("MLLW Tide Prediction (m)") +
  scale_x_datetime(date_breaks = "30 min", date_labels = "%H:%M", expand = c(0, 0)) +
  theme(plot.margin = margin(b = 0))

# Stitch plots together using patchwork
p_patch <- p1 / (p2 & theme(plot.margin = margin(0, 0, 0, 0, "pt"))) & ylab(NULL) & theme(plot.background = element_rect(fill = "transparent", colour = NA_character_))

# Use the tag label as a y-axis label
p_patch <- wrap_elements(p_patch) +
  labs(tag = "MLLW Tide Prediction (m)") +
  theme(
    plot.tag = element_text(size = rel(1.2), angle = 90),
    plot.tag.position = "left",
    plot.background = element_rect(fill = "transparent", colour = NA_character_)
  )


if (!dir.exists("Figures")) {
  dir.create("Figures")
}
path <- "Figures/Figure3."
save_pdf_png(plot = p_patch, path = path, w = 169, h = 90, u = "mm", bg = "transparent", dpi = 600)

First component of Figure 3 (second section added once data has been loaded

At three time points (see full Figure 3, finalized below) at the beginning and end of the sampling period as well as at low tide (at 14:00 PST), we collected triplicate 1 L samples from each location as biological replicates. At the beginning and end of the sampling period, we also filtered 1 L MilliQ water via the procedure described above to serve as negative field controls. Additionally, using an Orion Model 1230 meter (Orion Research Inc., Beverly, MA, USA), we recorded temperature and salinity in each location directly after samples were collected.

Once finished filtering, Sterivex cartridges were dried by pushing air through them using a sterile 3 mL syringe, capped, placed in sterile Whirl-Pak bags (Whirl-Pak, Madison, WI, USA). Then, samples were stored in a cooler on ice until transported back to the laboratory at the end of the sampling period. Samples were transferred to a -20°C freezer for up to 18 days, at which time they were processed to extract nucleic acids from the captured materials.

2.4 DNA Extraction & Library Preparation

Within 18 days of collection, we extracted DNA from the Sterivex cartridge using the DNeasy Blood and Tissue Kit (Qiagen, Germantown, MD, USA) and the modifications described in Spens et al. (2017). In short, we incubated the filter cartridge overnight with proteinase K and ATL. Then, we extracted the liquid from the cartridge with a syringe and mixed it with equal volumes of AL buffer and 0°C ethanol before proceeding with the manufacturer’s extraction protocol. One negative extraction control (DNA-grade water in place of a sample) was included in each of four batches of extractions. Extracted samples to be PCR amplified were stored at -20°C for up to 6 months.

We PCR amplified the extracted DNA in triplicate within 6 months of storage, using the mlCOIintF/ jgHCO2198 primer set targeting a 313 bp fragment of the mitochondrial COI region optimized by Leray et al. (2013) (Forward: GGWACWGGWTGAACWGTWTAYCCYCC; Reverse: TAIACYTCIGGRTGICCRAARAAYCA) with Nextera modifications. Following Curd et al. (2019), we used a 25 \(\mu\)l PCR reaction mixture consisting of 12.5 \(\mu\)l of Qiagen Multiplex Mix (Qiagen, Germantown, MD, USA), 2.5 \(\mu\)l each of the forward and reverse primers at 2 \(\mu\)M (Integrated DNA Technologies, Inc., Coralville, IA, USA), 6.5 \(\mu\)l of PCR-grade water, and 1 \(\mu\)l of undiluted DNA template. The PCR thermocycling touchdown profile began with an initial denaturation at 95° for 15 minutes to activate the DNA polymerase, followed by 13 cycles of denaturation (94° for 30 seconds), annealing (starting at 69.5° for 30 seconds, with the temperature decreased by 1.5° each cycle), and extension (72° for 1 minute). Then, an additional 35 cycles were run with the same denaturation and extension steps as above with an annealing temperature of 50°, followed by a final extension at 72° for 10 minutes. PCR reactions were prepared in a designated DNA-free hood until the template was added.

PCR amplification was conducted in two batches; in each batch, we included one no-template negative PCR control (DNA-grade water used as template). Additionally, we extracted DNA from tissue from five organisms across a range of phyla we expected to amplify with the mlCOIintF/ jgHCO2198 primers, but not expected to be present at Pillar Point in particular (Mytilus edulis, Mizuhopecten yessoensis, Xiphias gladius, Mercenaria mercenaria, Lutjanus campechanus) using the standard tissue extraction protocol detailed in the DNeasy Blood and Tissue Kit (Qiagen, Germantown, MD, USA). These tissues were obtained from a local grocery store and it was assumed that they were labeled correctly, although previous work has indicated mislabeling in seafood stores can occur (e.g. Willette et al., 2017). Extracts from the 5 tissue samples were combined in equimolar amounts to form a mock community used as a positive PCR control in each batch. Triplicate PCR amplicons from both samples and controls were not subsequently pooled, but were carried through the remaining library preparation and sequencing steps as technical replicates. We electrophoresed and visualized a subset of PCR products on a 1.5% agarose gel stained with GelRed® (Biotium, Fremont, CA, USA) to ensure successful amplification and correct product sizes as well as lack of contamination.

Post-PCR library preparation and sequencing was conducted at the Georgia Genomics and Bioinformatics Core (GGBC, UG Athens, GA, RRID:SCR_010994). In short, provided PCR amplicons were cleaned using AMPure XP magnetic beads (Beckman Coulter, Indianapolis, IN, USA), barcoded with Nextera adapters (Illumina, San Diego, CA, USA) during a second PCR (3 min at 95°C; 15 cycles of 30 sec at 95°C, 30 sec at 67°C, and 30 sec at 72°C; and 4 min at 72°C), cleaned again using AMPureXP magnetic beads, and pooled in equimolar ratio. The resulting library was sequenced on a MiSeq PE 2x250bp (500 cycles) using Reagent Kit V2 with 25% PhiX spike-in (Illumina, San Diego, CA, USA). Given our technical replication of samples and controls, our final library included 6 negative field controls (1 at beginning and end of field sampling, amplified in triplicate), 12 negative extraction controls (1 in each of 4 extraction sets, amplified in triplicate), 2 positive PCR controls (1 in each of 2 amplification batches), 2 no-template negative PCR controls (1 in each of 2 amplification batches), and 159 samples (53 field samples, amplified in triplicate).

2.5 Bioinformatics

We processed sequencing data using the Anacapa Toolkit, which contains two core modules: one for quality control and ASV parsing, and one for classifying taxonomy (Curd et al., 2019). Briefly, we ran the first module using default parameters, which uses cutadapt (version 1.16) (Martin, 2011) for adapter and primer trimming, FastX-Toolkit (version: 0.0.13) (citation) for quality trimming, and dada2 (version 1.6) (Callahan et al., 2016) for assigning ASVs. For the second module, we utilized the MIDORI2 reference database, a quality controlled and updated database built from GenBank release 253 (20 December 2022) that has been technically validated (Leray et al., 2022). Following Gold et al. (2022), we adjusted the identity and query coverage to 95% (default: 80%) to account for the relative incompleteness of the broad COI reference database compared to more taxonomically-specific databases (Curd et al., 2019). The second module relies on Bowtie 2 (version 2.3.5) (Langmead and Salzberg, 2012) and a modified instance of BLCA (Gao et al., 2017) as dependencies. Following Gold et al. (2021a) we only kept taxonomic assignments that had a bootstrap confidence cutoff score of 60 or higher in BLCA, to avoid spurious assignments from the incomplete reference database. We modified the Anacapa Container (Ogden, 2018), a Singularity container with all the needed dependencies for executing the Anacapa Toolkit, to enable the pipeline to be run in a high-performance computing environment requiring two-step authentication; the updated container, scripts, and reference database with the required Bowtie 2 index library needed to reproduce our bioinformatics process can be found on Dryad (link/citation).

Raw ASVs, taxonomy assignments, and metadata were converted into interchangeable ampvis2 (version 2.8.6 packageVersion("ampvis2")) (Andersen et al., 2018) and phyloseq (version 1.46.0 packageVersion("phyloseq")) (McMurdie and Holmes, 2013) objects in R version 4.3.1 (2023-06-16) R.Version()$version.string to facilitate decontamination and further analyses. Singletons were removed using ampvis2, and samples were further decontaminated using phyloseq by removing all ASVs that appeared in any negative field control, negative extraction control, or no-template negative PCR control, a choice made due to the very low number of overlapping ASVs between samples and negative controls (Table S1). Samples were rarified to the minimum number of reads of any sample using ampvis2. Subsequent analyses were robust to all three decontamination steps.

### Importing Anacapa###
# Read in Anacapa Pipeline Output and process into ampvis2 objects for use in analyses

# Read in Anacapa Output
# Change p_confidence to switch to a different percent confidence output
###
p_confidence <- 60
###
current_subfolder <- "Data/Anacapa Pipeline Output/Tide Pool 5184 - 31 January 2023"

filepath <- paste0(
  current_subfolder,
  "/CO1/CO1_taxonomy_tables/Summary_by_percent_confidence/",
  p_confidence,
  "/CO1_ASV_raw_taxonomy_",
  p_confidence,
  ".txt"
)

AnacapaOutput <- read.table(filepath, header = TRUE, sep = "\t", dec = ".")

# Format ASV Table (extract just the needed information)
endcol <- ncol(AnacapaOutput) - 1
asvmat <- AnacapaOutput[, 2:endcol]
rownames(asvmat) <- AnacapaOutput[, 1]

# Format Taxonomy Table (extract just the needed information and rename columns)
endcol <- ncol(AnacapaOutput)
taxmat <- AnacapaOutput[, endcol]
taxmat <- str_split_fixed(taxmat, ";", 7)
colnames(taxmat) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
rownames(taxmat) <- rownames(asvmat)

# Read Metadata (saved in multiple different files)

current_subfolder <- paste0(current_subfolder, "/Metadata-Added/PillarPoint_SampleInfo_")

filepath <- paste0(current_subfolder, "EnviroData.txt")
EnviroData <- read.table(filepath, header = TRUE, sep = "\t", dec = ".", strip.white = TRUE)

filepath <- paste0(current_subfolder, "FieldSampling.txt")
FieldSampling <- read.table(filepath, header = TRUE, sep = "\t", dec = ".", strip.white = TRUE)
FieldSampling <- FieldSampling[, 1:4]

filepath <- paste0(current_subfolder, "Sequencing.txt")
Sequencing <- read.table(filepath, header = TRUE, sep = "\t", dec = ".", strip.white = TRUE)

filepath <- paste0(current_subfolder, "Location.txt")
Location <- read.table(filepath, header = TRUE, sep = "\t", dec = ".", strip.white = TRUE)

# Format Metadata
sampledata <- merge(FieldSampling, EnviroData, all = T)
sampledata <- sampledata %>% drop_na(Dereplicated_Sample_Name)
# removing one sample that has EnviroData reading but the field sample was not processed

sampledata <- merge(Sequencing, sampledata, all = T)
rownames(sampledata) <- sampledata[, "X"]

sampledata$Dereplicated_Sample_Name <- as.factor(sampledata$Dereplicated_Sample_Name)
sampledata$Type <- as.factor(sampledata$Type)
sampledata$PCR <- as.factor(sampledata$PCR)
sampledata$Time <- as.factor(sampledata$Time)
sampledata$Location <- as.factor(sampledata$Location)
sampledata$Extraction <- as.factor(sampledata$Extraction)
sampledata$SiteByTime <- as.factor(paste(sampledata$Location, sampledata$Time))

# Make ampvis2 objects without any pre-processing
asvmat_amp <- tibble::rownames_to_column(asvmat, "ASV")
sampledata_amp <- tibble::rownames_to_column(sampledata, "Sample")
taxmat_amp <- tibble::rownames_to_column(as.data.frame(taxmat), "ASV")
amp <- amp_load(asvmat_amp, metadata = sampledata_amp, taxonomy = taxmat_amp)

amp_nocontrols <- amp_filter_samples(amp, Type %in% c("sample"))

## 22 samples and 3315 OTUs have been filtered 
## Before: 181 samples and 71947 OTUs
## After: 159 samples and 68632 OTUs

amp_controls <- amp_filter_samples(amp, Type %in% c("control"))

## 159 samples and 68607 OTUs have been filtered 
## Before: 181 samples and 71947 OTUs
## After: 22 samples and 3340 OTUs

# Make ampvis2 objects with singletons removed
amp_NS <- amp_load(asvmat_amp, metadata = sampledata_amp, taxonomy = taxmat_amp, pruneSingletons = TRUE)

## 41532 singletons have been removed

amp_nocontrols_NS <- amp_filter_samples(amp_NS, Type %in% c("sample"))

## 22 samples and 381 OTUs have been filtered 
## Before: 181 samples and 30415 OTUs
## After: 159 samples and 30034 OTUs

amp_controls_NS <- amp_filter_samples(amp_NS, Type %in% c("control"))

## 159 samples and 30009 OTUs have been filtered 
## Before: 181 samples and 30415 OTUs
## After: 22 samples and 406 OTUs

# Make ampvis2 objects with a full decontamination procedure

# Function that removes any ASV that appears in any control from all samples
full_decontamination <- function(amp) {
  physeq <- amp_to_phyloseq(amp, NULL, NULL, NULL)

  physeq_controls <- physeq %>%
    subset_samples(Type == "control") %>%
    prune_taxa(taxa_sums(.) > 0, .)

  badASV <- taxa_names(physeq_controls)
  allASV <- taxa_names(physeq)
  allASV <- allASV[!(allASV %in% badASV)]
  physeq_FD <- prune_taxa(allASV, physeq)

  av2_otutable <- data.frame(otu_table(physeq_FD)@.Data)
  av2_taxtable <- data.frame(tax_table(physeq_FD)@.Data)
  av2_metadata <- data.frame(sample_data(physeq_FD))

  amp_return <- amp_load(av2_otutable, metadata = av2_metadata, taxonomy = av2_taxtable)
  amp_return <- amp_filter_samples(amp_return, Type %in% c("sample"))

  return(amp_return)
}

amp_nocontrols_FD <- full_decontamination(amp)

## 22 samples and 0 OTUs have been filtered 
## Before: 181 samples and 68607 OTUs
## After: 159 samples and 68607 OTUs

amp_nocontrols_NS_FD <- full_decontamination(amp_NS)

## 22 samples and 0 OTUs have been filtered 
## Before: 181 samples and 30009 OTUs
## After: 159 samples and 30009 OTUs

# Make ampvis2 objects by rarefying to the minimum number of reads of any sample

# Function that wraps around an ampvis2 function to rarefy samples to the minimum number of reads of any sample
rarefaction <- function(amp) {
  set.seed(0)
  reads <- colSums(amp$abund)
  minreads <- min(reads)
  amp_R <- amp_subset_samples(amp, rarefy = minreads)
  return(amp_R)
}

amp_nocontrols_R <- rarefaction(amp_nocontrols)

## 0 samples have been filtered.

amp_nocontrols_NS_R <- rarefaction(amp_nocontrols_NS)

## 0 samples have been filtered.

amp_nocontrols_FD_R <- rarefaction(amp_nocontrols_FD)

## 0 samples have been filtered.

amp_nocontrols_NS_FD_R <- rarefaction(amp_nocontrols_NS_FD)

## 0 samples have been filtered.

### Making GBIF Data Frames###
## Create data frames for submission to GBIF
# Export version of data with decontamination and controls removed but no additional processing

GBIF_export <- amp_nocontrols_FD

GBIF_export$metadata <- merge(GBIF_export$metadata, Location, all = T)

TimeZone <- "UTC−08:00"
Date <- "28 January 2022"
GBIF_export$metadata <- cbind(GBIF_export$metadata, TimeZone)
GBIF_export$metadata <- cbind(GBIF_export$metadata, Date)

names(GBIF_export$metadata)[names(GBIF_export$metadata) == "Sample"] <- "Sequencing_Sample_Name"

col_order <- c("Sequencing_Sample_Name", "Sample_Name", "Location", "Latitude", "Longitude", "CRS", "Date", "Time", "TimeZone", "T", "S")
GBIF_export$metadata <- GBIF_export$metadata[, col_order]

# Add sequences to file

current_subfolder <- "Data/Anacapa Pipeline Output/Tide Pool 5184 - 31 January 2023"
filepath <- paste0(
  current_subfolder,
  "/CO1/CO1_taxonomy_tables/CO1_ASV_taxonomy_detailed.txt"
)

AnacapaDetailed <- read.table(filepath, header = TRUE, sep = "\t", dec = ".")
sequences <- as.data.frame(AnacapaDetailed[, c("sequence", "sequencesF", "sequencesR")])
rownames(sequences) <- AnacapaDetailed$CO1_seq_number

GBIF_export$tax <- merge(sequences, GBIF_export$tax, by = "row.names", all.x = FALSE, all.y = TRUE, sort = FALSE)
rownames(GBIF_export$tax) <- GBIF_export$tax$Row.names

col_order <- c("sequence", "sequencesF", "sequencesR", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
GBIF_export$tax <- GBIF_export$tax[, col_order]



if (!dir.exists("Analysis Products/Processed eDNA/GBIF")) {
  dir.create("Analysis Products/Processed eDNA/GBIF", recursive = TRUE)
}
write.csv(GBIF_export$abund, "Analysis Products/Processed eDNA/GBIF/ASVTable.csv", row.names = TRUE)
write.csv(GBIF_export$tax, "Analysis Products/Processed eDNA/GBIF/TaxTable.csv", row.names = TRUE)
write.csv(GBIF_export$metadata, "Analysis Products/Processed eDNA/GBIF/SampleData.csv", row.names = FALSE)

To ensure the accuracy of taxonomic assignments, we analyzed occurrence data from the Global Biodiversity Information Facility to investigate whether identified species had occurrence records in the California Current System and known ranges that encompassed Pillar Point, using the spocc (version 1.2.2 packageVersion("spocc")) (Chamberlain, 2021) package. A phylogenetic tree based on taxonomic assignments was created using the taxize (version 0.9.100 packageVersion("taxize")) (Chamberlain and Szöcs, 2013) and ggtreeExtra 1.12.0 packageVersion("ggtreeExtra") (citation) packages.

Data Analysis

To understand whether eDNA signals could be distinguished by location, we first analyzed individual-level differences between locations; that is, whether ASVs or individual taxa (agglomerated species-level taxonomic assignments) were unique to, or associated with, particular locations. We calculated and visualized unique ASVs and taxa using the eulerr (version 7.0.0 packageVersion("eulerr")) (Larsson, 2022) package. However, with metabarcoding data in particular, taxa or ASVs that are unique to a given location are not necessarily ecologically meaningful; they could include rare taxa present elsewhere but not amplified and exclude taxa that are well correlated with particular locations but sometimes detected at others. Thus, we also analyzed ASVs and taxa using an indicator species framework (Dufrêne and Legendre, 1997). In this framework, the null hypothesis is that the frequency of a taxon’s or ASV’s presence in samples from a particular location is not higher that the frequency of that taxon’s or ASV’s presence in samples from other locations. For each location, we identified all statistically-significant indicator taxa and ASVs with indicator value indices above 0.7—that is, indicators that are well-associated with a site group even if they are detected in samples from other sites—based on presence-absence data per sample using the indicspecies package (version 1.7.14 packageVersion("indicspecies")) (Cáceres and Legendre, 2009).

We followed our individual-level analyses with approaches to test whether community composition varied by location. We calculated a Jaccard dissimilarity matrix across all samples using vegan (version 2.6.4 packageVersion("vegan")) (Oksanen et al., 2022). Then, we tested for differences in community composition among locations using a permutational multivariate analysis of variance (PERMANOVA) with the model eDNA Presence ~ Location + Time + Biological Replicates using the adonis function in vegan. We confirmed that locational differences found via PERMANOVA were not a result of differences in dispersions by testing for homogeneity of dispersions using the betadisper function in vegan. We also visualized Jaccard dissimilarity using non-metric multidimensional scaling (NMDS) using the metaMDS function in vegan. We coupled these analyses with a partitioning among medoids algorithm (Kaufman and Rousseeuw, 1990) implemented with the pam function in the cluster package (version 2.1.6 packageVersion("cluster")) (Maechler et al., 2022). Rather than assuming location clusters a priori, we validated the optimal number of clusters for a given dataset by finding the number of clusters (k) that maximized the average silhouette width, a measure of how well-structured the clusters are. Finally, to better understand how the difference between locations varied over the time sampled, for each time point, we calculated the Jaccard dissimilarity between each unique combination of replicates within each site, and then the pairwise Jaccard dissimilarity between each unique combination of samples across the three pairs of sites: S1-N, S2-N, and S1-S2.

To further understand whether differences in eDNA detections corresponded with underlying ecological gradients, we compared the unique and indicator taxa identified at each site to their ecological zonation in a highly regarded field guide to the Pacific intertidal, Between Pacific Tides (Ricketts et al., 1985). To account for potential variations in taxonomic names between Between Pacific Tides and the MIDORI2 reference database, we used the World Register of Marine Species (WoRMS) to identify all synonymized names for each unique and indicator taxon identified, and we searched Between Pacific Tides for all synonyms. We compared the proportion of high and middle intertidal species identified across each location using a chi-squared test with p-values computed via Monte Carlo simulation.

Results

3.1 Sequencing Results & Taxonomic Diversity

### Full Output Metrics###
# Generate information about the full Anacapa Pipeline results

# Gather information about full Anacapa Pipeline results before any processing
total_ASVs <- nrow(amp$abund)
total_reads <- sum(amp$abund)

samples <- length(amp$metadata$Type[amp$metadata$Type == "sample"])
controls <- length(amp$metadata$Type[amp$metadata$Type == "control"])

# Check positive PCR controls
amp_MC <- amp_filter_samples(amp, Location %in% c("MC"))

## 179 samples and 71797 OTUs have been filtered 
## Before: 181 samples and 71947 OTUs
## After: 2 samples and 150 OTUs

df <- merge(amp_MC$abund, amp_MC$tax, by = "row.names")
df <- df[, c(2:3, 10)]
colnames(df) <- c("MC1", "MC2", "Species")
agg_df <- df %>%
  group_by(Species) %>%
  summarise(MC1_Reads = sum(MC1), MC2_Reads = sum(MC2)) %>%
  arrange(desc(MC1_Reads))
agg_df[agg_df == ""] <- "No Taxonomic Assignment"
agg_df$Species <- remove_NCBInums(agg_df$Species)
agg_df

amp_no_MC <- amp_filter_samples(amp, !(Location %in% c("MC")))

## 2 samples and 150 OTUs have been filtered 
## Before: 181 samples and 71947 OTUs
## After: 179 samples and 71797 OTUs

if (length(intersect(rownames(amp_MC$abund), rownames(amp_no_MC$abund))) == 0) {
  MC_overlap <- "no ASVs present in the positive PCR controls occurred in any other samples, providing no evidence of index hopping"
} else {
  MC_overlap <- "some ASVs present in the positive PCR controls occurred in other samples, providing evidence of index hopping (CHECK THIS)"
}

# Check intersection between samples and controls
intersect <- intersect(rownames(amp_controls$abund), rownames(amp_nocontrols$abund))
overlap <- length(intersect)

overlap_table <- data.frame(ASV = intersect, "Taxonomic Assignment" = amp_nocontrols$tax[amp_nocontrols$tax$OTU %in% intersect, ]$Species)
overlap_table$Taxonomic.Assignment <- remove_NCBInums(overlap_table$Taxonomic.Assignment)
overlap_table

Using the Anacapa Toolkit, we identified 71947 (total_ASVs) ASVs from 9291419(total_reads) reads across 159 (samples) samples and 22 (controls) controls (positive PCR controls, negative field controls, negative extraction controls, and no-template negative PCR controls). All expected taxa amplified in each positive PCR control, and no ASVs present in the positive PCR controls occurred in any other samples, providing no evidence of index hopping (MC_overlap). Before application of any decontamination steps, only 25 (overlap) ASVs were shared between the samples and negative field, extraction, and PCR controls.

### Core Data Sets###
# Establish processed data set used throughout analyses

data_amp <- amp_nocontrols_NS_FD_R

data_physeq <- amp_to_phyloseq(data_amp, NULL, NULL, NULL)

data_amp_taxa <- data_amp
data_amp_taxa$abund <- aggregate_abund(data_amp$abund, data_amp$tax, tax_aggregate = "Species", format = "abund")

## 24191 OTUs (out of 27038) with no assigned taxonomy at Species level were removed before aggregating OTUs

if (!dir.exists("Analysis Products/Processed eDNA/Manuscript Analysis")) {
  dir.create("Analysis Products/Processed eDNA/Manuscript Analysis", recursive = TRUE)
}
write.csv(data_amp$abund, "Analysis Products/Processed eDNA/Manuscript Analysis/ASVTable.csv", row.names = TRUE)
write.csv(data_amp$tax, "Analysis Products/Processed eDNA/Manuscript Analysis/TaxTable.csv", row.names = TRUE)
write.csv(data_amp$metadata, "Analysis Products/Processed eDNA/Manuscript Analysis/SampleData.csv", row.names = FALSE)

### Figure 3--Continued###
# Finish the bottom portion of Figure 3 now that data has been loaded

data_forchart <- amp$metadata[!is.na(amp$metadata$Time), ]

data_forchart <- data_forchart %>%
  group_by(SiteByTime) %>%
  mutate(chart_num = ceiling(row_number() / 3))

facet.labs <- c("Tide Pool 1\n (S1)", "Tide Pool 2\n (S2)", "Nearshore\n (N)", "Field Blank")
names(facet.labs) <- c("S1", "S2", "N", "FB")

data_forchart$Location <- factor(data_forchart$Location, levels = c("S1", "S2", "N", "FB"))
data_forchart$chart_num <- factor(data_forchart$chart_num, levels = c("3", "2", "1"))

sample_schema <- ggplot(data = data_forchart, aes(
  x = Time, y = chart_num, shape = Location,
  color = Time,
  fill = after_scale(alpha(colour, 0.05))
)) +
  geom_beeswarm(cex = 2, size = 3) +
  geom_box(
    aes(
      xmin = after_stat(x) - 0.45,
      xmax = after_stat(x) + 0.45,
      ymin = after_stat(y) - 0.45,
      ymax = after_stat(y) + 0.45
    ),
    radius = unit(5, "pt")
  ) +
  scale_shape_manual(values = c(15, 17, 16, 18)) +
  scale_color_viridis_d(direction = -1, begin = 0, end = .97) +
  facet_grid(Location ~ .,
    scales = "free_y", space = "free_y",
    switch = "y", labeller = labeller(Location = facet.labs)
  ) +
  theme_void(base_size = 12) +
  theme(
    legend.position = "none",
    strip.text.y.left = element_text(angle = 0, hjust = 1, vjust = .9),
    plot.margin = margin(0, 0, 0, 0, "pt")
  )


combined <- (plot_spacer() + p_patch + plot_spacer() + plot_layout(widths = c(3, 50, 1))) / (plot_spacer() + sample_schema + plot_layout(widths = c(1, 50))) & theme(plot.margin = margin(0, 0, 0, 0, "pt"))


path <- "Figures/Figure3_full."
save_pdf_png(plot = combined, path = path, w = 169, h = 169, u = "mm", bg = "transparent", dpi = 600)

Figure 3: Graph of tide predictions during field sampling on 28 January 2022, with icons to indicate the sampling scheme

### Processed Output Metrics###
# Gather information about processed Anacapa Toolkit

processed_ASVs <- nrow(data_amp$abund)
processed_reads <- sum(data_amp$abund)

# Create a phyloseq object where all samples are combined into one
total <- as.data.frame(rowSums(data_amp$abund))
physeq_summed <- amp_to_phyloseq(NULL, total, data_amp$tax[, 1:7], NULL)

# Summarize the number of unassigned ASVs
kingdom <- as.data.frame(tax_table(physeq_summed)[, 1])
unassigned_ASVs <- kingdom %>%
  group_by(Kingdom) %>%
  summarise(n = n()) %>%
  mutate(freq = n / sum(n))

# Summarize number of unique taxa at each taxonomic level
unique <- apply(tax_table(physeq_summed), 2, function(x) length(unique(x[x != ""])))

### GBIF Set-Up###
# Create functions and datasets needed for GBIF analysis

# Create list of identified species
species_list <- unique(data_amp$tax$Species)
species_list <- remove_NCBInums(species_list)
species_list <- species_list[species_list != ""]
species_list <- sub("_", "-", species_list) # to address Pseudo_nitzchia error

# Standardize species names against GBIF backbone
species_GBIF <- name_backbone_checklist(species_list)
species_standardized <- species_GBIF$species[!is.na(species_GBIF$species)]
species_both <- species_GBIF[c("species", "verbatim_name")]

## Core GBIF checking function
distribution_check <- function(small_bound, big_bound, species, min_records, efficiency) {
  ## Function that checks GBIF given a particular geographic region
  boundary_function <- function(species, bounds) {
    no_data <- c()
    confirmed <- c()
    for (i in species) {
      print(i)
      occ <- occ(
        query = i,
        geometry = bounds,
        from = "gbif",
        limit = min_records,
        has_coords = TRUE,
        gbifopts = list(hasGeospatialIssue = FALSE)
      )
      occ.df <- occ2df(occ)
      if (nrow(occ.df) >= min_records) {
        confirmed <- append(confirmed, i)
      } else {
        no_data <- append(no_data, i)
      }
    }
    return(list(no_data = no_data, confirmed = confirmed))
  }

  ## Function that checks GBIF to confirm whether a particular location is within the range of observations
  range_function <- function(species) {
    CC.confirmed <- c()
    CC.not_confirmed <- c()
    CC.no_data <- c()
    limits <- c(10, 100, 1000, 10000)

    for (i in species) {
      track <- 0
      for (l in limits) {
        print(i)
        print(l)
        occ <- occ(i, from = "gbif", limit = l, has_coords = TRUE, gbifopts = list(hasGeospatialIssue = FALSE))
        occ.df <- occ2df(occ)
        if (nrow(occ.df) < min_records) {
          CC.no_data <- append(CC.no_data, i)
          print("no data")
          track <- 1
          break
        }
        bb_occ <- sp::bbox(cbind(occ.df$longitude, occ.df$latitude))
        if (bb_occ["x", "min"] < small_bound[1] &
          bb_occ["x", "max"] > small_bound[1] &
          bb_occ["y", "min"] < small_bound[4] &
          bb_occ["y", "max"] > small_bound[4]) {
          CC.confirmed <- append(CC.confirmed, i)
          print("confirmed")
          track <- 1
          break
        }
        if (nrow(occ.df) < l) {
          CC.not_confirmed <- append(CC.not_confirmed, i)
          print("not confirmed")
          track <- 1
          break
        }
      }
      if (track == 0) {
        CC.not_confirmed <- append(CC.not_confirmed, i)
        print("not confirmed")
      }
    }
    return(list(no_data = CC.no_data, not_confirmed = CC.not_confirmed, confirmed = CC.confirmed))
  }


  # Actual Steps
  PP <- boundary_function(species, small_bound)

  if (efficiency == TRUE) {
    CC <- boundary_function(PP[["no_data"]], big_bound)
  } else {
    CC <- boundary_function(species, big_bound)
  }

  if (efficiency == TRUE) {
    range_confirmed <- range_function(CC[["confirmed"]])
    range_unconfirmed <- range_function(CC[["not_confirmed"]])
    return(list(PP, CC, range_confirmed, range_unconfirmed))
  } else {
    range <- range_function(species)
    return(list(PP, CC, range))
  }
}

# Get the California Current System shape

CC_shape <- mr_features_get(type = "MarineRegions:lme", featureID = "lme.10", format = "json")

# check that it's the right polygon
class(CC_shape) <- "mr_geojson"
CC_shape <- mr_as_wkt(CC_shape)
pnt <- st_as_sfc(CC_shape, crs = 4326)
mapview(pnt)

# Simplify the polygon for use with GBIF
CC_shape_sp <- st_transform(pnt, crs = "+proj=aeqd +lat_0=53.6 +lon_0=12.7")
mapview(CC_shape_sp)

CC_shape_sp <- st_simplify(CC_shape_sp, preserveTopology = FALSE, dTolerance = 10000)
CC_shape_sp <- st_transform(CC_shape_sp, crs = 4326)
mapview(CC_shape_sp)

### GBIF Run###
# Run function that does GBIF analysis (this takes a long time; recommend running with job::job() when working in chunks)

distribution <- distribution_check(c(-122.50, 37.49, -122.49, 37.50), CC_shape_sp, species_standardized, 5, FALSE)

### GBIF Summary###
# Organize the GBIF analysis output

d <- distribution
PP <- d[[1]]$confirmed
CC <- d[[2]]$confirmed
range <- d[[3]]$confirmed
no_data <- d[[3]]$no_data

df_PP <- c()
df_CC <- c()
df_range <- c()
df_no_data <- c()
df_not_in_GBIF <- c()


for (s in species_both$species) {
  if (is.na(s)) {
    df_not_in_GBIF <- append(df_not_in_GBIF, TRUE)
    df_PP <- append(df_PP, NA)
    df_CC <- append(df_CC, NA)
    df_range <- append(df_range, NA)
    df_no_data <- append(df_no_data, NA)
    next
  } else {
    df_not_in_GBIF <- append(df_not_in_GBIF, FALSE)
  }

  if (s %in% PP) {
    df_PP <- append(df_PP, TRUE)
  } else {
    df_PP <- append(df_PP, FALSE)
  }

  if (s %in% CC) {
    df_CC <- append(df_CC, TRUE)
  } else {
    df_CC <- append(df_CC, FALSE)
  }

  if (s %in% range) {
    df_range <- append(df_range, TRUE)
  } else {
    df_range <- append(df_range, FALSE)
  }

  if (s %in% no_data) {
    df_no_data <- append(df_no_data, TRUE)
  } else {
    df_no_data <- append(df_no_data, FALSE)
  }
}

GBIF_summary <- data.frame(
  Species = species_both$verbatim_name,
  GBIFSpecies = species_both$species,
  PillarPointOccurrence = df_PP,
  CCSOccurrence = df_CC,
  Range = df_range,
  NoData = df_no_data,
  NotInGBIF = df_not_in_GBIF
)

percentages <- GBIF_summary %>%
  group_by(across(PillarPointOccurrence:NotInGBIF)) %>%
  summarise(count = n()) %>%
  ungroup() %>%
  mutate(per = count / sum(count))

known_summary <- as.numeric(GBIF_summary %>% filter(!is.na(GBIFSpecies)) %>% summarise(count = sum(Range == TRUE & CCSOccurrence == TRUE)))

Focused on just the 159 (samples) samples, we removed singletons and all 25 (overlap) ASVs that overlapped with negative controls (Supplemental Information X), and we rarefied samples to the minimum number of reads of any sample, leaving 27038 (processed_ASVs) ASVs and 4047981 (processed_reads) reads. 89.1% (sprintf("%0.1f%%", unassigned_ASVs[unassigned_ASVs$Kingdom == "",]$freq * 100)) of ASVs were unassigned, and the remainder included representatives from 28 (unique["Phylum"]) phyla, 59 (unique["Class"]) classes, 132 (unique["Order"]) orders, 234 (unique["Family"]) families, 325 (unique["Genus"]) genera, and 415 (unique["Species"]) species (Figure 4). A breakdown of the relative proportions of reads and ASVs assigned to each phyla can be found in Table 1. To ensure the accuracy of these taxonomic identifications, we confirmed that 293 (known_summary) of 415 (unique["Species"]) identified species (70.6% sprintf("%0.1f%%", known_summary/unique["Species"] * 100)) have occurrence records in the California Current System (CCS) and known ranges that encompass Pillar Point. (additional details in SI 1 and Table S2).

### TABLE 1###
# Create table with summary of reads, ASVs, and species by phyla

# Remove all unclassified ASVs and condense to phylum level
physeq_totalreads <- tax_glom(physeq_summed, taxrank = "Phylum", NArm = TRUE)
Phyla_Reads <- data.frame(
  Phyla = tax_table(physeq_totalreads)[, 2],
  Reads = rowSums(otu_table(physeq_totalreads))
)

physeq_totalASVs <- physeq_summed

otu_table(physeq_totalASVs)[otu_table(physeq_totalASVs) > 0] <- 1
physeq_totalASVs <- tax_glom(physeq_totalASVs, taxrank = "Phylum", NArm = TRUE)


Phyla_ASVs <- data.frame(Phyla = tax_table(physeq_totalASVs)[, 2], ASVs = rowSums(otu_table(physeq_totalASVs)))


physeq_speciesbyphyla <- physeq_summed
physeq_speciesbytaxa <- tax_glom(physeq_summed, taxrank = "Species", NArm = TRUE)
otu_table(physeq_speciesbytaxa)[otu_table(physeq_speciesbytaxa) > 0] <- 1
physeq_speciesbytaxa <- tax_glom(physeq_speciesbytaxa, taxrank = "Phylum", NArm = TRUE)

Phyla_Species <- data.frame(Phyla = tax_table(physeq_speciesbytaxa)[, 2], Species = rowSums(otu_table(physeq_speciesbytaxa)))


list <- list(Phyla_Reads, Phyla_ASVs, Phyla_Species)
full <- list %>% reduce(full_join, by = "Phylum")
full[is.na(full)] <- 0


full_sorted <- full[order(-full$Reads), ]

full_sorted <- full_sorted %>%
  mutate(Reads_Percent = scales::label_percent()(Reads / sum(Reads))) %>%
  mutate(ASVs_Percent = scales::label_percent()(ASVs / sum(ASVs))) %>%
  mutate(Species_Percent = scales::label_percent()(Species / sum(Species)))

full_sorted$Phylum <- remove_NCBInums(full_sorted$Phylum)
full_sorted$Phylum <- str_replace(full_sorted$Phylum, "phylum_", "Class: ")

col_order <- c("Phylum", "Reads", "Reads_Percent", "ASVs", "ASVs_Percent", "Species", "Species_Percent")

full_sorted <- full_sorted[, col_order]


write.csv(full_sorted, "Analysis Products/Table1.csv", row.names = FALSE)

full_sorted

### Figure 4###
# Create phylogenetic tree for all taxonomic assignments

options(getClass.msg = FALSE)

# Aggregate to species assignment and remove the NCBI identifiers
pt_physeq <- tax_glom(data_physeq, taxrank = "Species", NArm = FALSE)
tax_table(pt_physeq) <- tax_table(pt_physeq)[!tax_table(pt_physeq)[, "Kingdom"] == ""]


tax_table(pt_physeq)[tax_table(pt_physeq) == ""] <- NA
tax_table(pt_physeq) <- tax_table(pt_physeq)[!is.na(tax_table(pt_physeq)[, "Phylum"])]
numbers <- sub(".*_", "", tax_table(pt_physeq))

taxa <- apply(numbers, 1, function(x) tail(na.omit(x), 1))

tax_table(pt_physeq) <- remove_NCBInums(tax_table(pt_physeq))


taxa_classification <- classification(taxa, db = "ncbi")


taxa_tree <- class2tree(taxa_classification, check = TRUE)

# Change the tip labels to match phyloseq objects (code only works if order is the same)
taxa_tree$phylo$tip.label <- row.names(tax_table(pt_physeq))


phyla <- unique(tax_table(pt_physeq)[, "Phylum"])

phyla_nodes <- c()
phyla_nodes <- data.frame("node_phylum" = c(), "node" = c())
phyla_nodes_secondary <- data.frame("node_phylum" = c(), "tips" = c())

for (x in phyla) {
  OTUs <- na.omit(row.names(tax_table(pt_physeq))[tax_table(pt_physeq)[, "Phylum"] == x])
  OTUs <- as.vector(OTUs)
  num <- ape::getMRCA(taxa_tree$phylo, tip = OTUs)
  if (is.null(num)) {
    nodes <- as_tibble(taxa_tree$phylo)
    num <- nodes$node[nodes$label == OTUs]
  }
  df <- data.frame("node_phylum" = c(x), "node" = c(num))
  phyla_nodes <- rbind(phyla_nodes, df)

  if (length(OTUs) < 6) {
    for (y in OTUs) {
      df2 <- data.frame("node_phylum" = c(x), "tips" = c(y))
      phyla_nodes_secondary <- rbind(phyla_nodes_secondary, df2)
    }
  }
}

pt_physeq <- merge_phyloseq(pt_physeq, taxa_tree$phylo)

# Use Polychrome package to create a good color palette for the number of phyla
set.seed(567629)
colors <- createPalette(28, c("#FFFFFF"), range = c(30, 80), target = c("normal", "protanope", "deuteranope", "tritanope"))

# Check whether the given order of colors has good contrast for different forms of colorblindness
pal.safe(colors)

# Manually re-order the colors to ensure there is good contrast between adjacent colors
colors_ordered <- colors[c(
  "NC2",
  "NC11",
  "NC27",
  "NC25",
  "NC4",
  "NC5",
  "NC8",
  "NC24",
  "NC9",
  "NC10",
  "NC12",
  "NC13",
  "NC1",
  "NC14",
  "NC15",
  "NC16",
  "NC7",
  "NC17",
  "NC18",
  "NC19",
  "NC20",
  "NC21",
  "NC22",
  "NC23",
  "NC26",
  "NC6",
  "NC3",
  "NC28"
)]

# Confirm that the new order works well
pal.safe(colors_ordered)

# Create ggtree
p <- ggtree::ggtree(pt_physeq, layout = "circular")

ordered_names <- ggtree::get_taxa_name(p)
ordered_names <- unique(ordered_names)
ordered_names[match(rownames(tax_table(pt_physeq)), ordered_names)] <- tax_table(pt_physeq)[, "Phylum"]
ordered_names <- unique(ordered_names)

names(colors_ordered) <- ordered_names

p <- p + geom_fruit(geom = geom_tile, mapping = aes(fill = Phylum), width = 3, show.legend = FALSE) +
  scale_fill_manual(values = colors_ordered) +
  scale_color_manual(values = colors_ordered)

p <- p + geom_fruit(data = phyla_nodes_secondary, geom = geom_tile, mapping = aes(y = tips, fill = node_phylum), width = 6, show.legend = FALSE)
p <- p + ggtree::geom_hilight(data = phyla_nodes, mapping = aes(node = node, fill = node_phylum), show.legend = FALSE)

p <- p + ggtree::geom_cladelab(data = phyla_nodes, mapping = aes(node = node, label = node_phylum), geom = "text", angle = "auto", offset = 5, fontsize = 2, barcolor = NA, show.legend = FALSE)

p <- p +
  theme(
    panel.background = element_rect(fill = "transparent"),
    plot.background = element_rect(fill = "transparent", color = NA),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.background = element_rect(fill = "transparent"),
    legend.box.background = element_rect(fill = "transparent")
  )

path <- "Figures/Figure4."
save_pdf_png(plot = p, path = path, w = 169, h = 169, u = "mm", bg = "transparent", dpi = 600)

Figure 4: : Phylogenetic tree colored by the phyla identified across all eDNA samples. Italicized taxonomic names are those with no phylum-level assignment, so class-level assignments are used to differentiate (additional information added after R figure production)

3.2 Individual-Level Differences Between Locations

### Euler Diagrams###

eulerr_options(quantities = list(fontsize = 10), legend = list(fontsize = 14))

# Subset ampvis2 data by location

amp_S1 <- amp_filter_samples(data_amp, Location %in% c("S1"))

## 108 samples and 14365 OTUs have been filtered 
## Before: 159 samples and 27038 OTUs
## After: 51 samples and 12673 OTUs

amp_S2 <- amp_filter_samples(data_amp, Location %in% c("S2"))

## 105 samples and 16586 OTUs have been filtered 
## Before: 159 samples and 27038 OTUs
## After: 54 samples and 10452 OTUs

amp_N <- amp_filter_samples(data_amp, Location %in% c("N"))

## 105 samples and 15230 OTUs have been filtered 
## Before: 159 samples and 27038 OTUs
## After: 54 samples and 11808 OTUs

# List the unique taxa at each location

S1_taxa <- unique(amp_S1$tax$Species)
S1_taxa <- S1_taxa[S1_taxa != ""]
S2_taxa <- unique(amp_S2$tax$Species)
S2_taxa <- S2_taxa[S2_taxa != ""]
N_taxa <- unique(amp_N$tax$Species)
N_taxa <- N_taxa[N_taxa != ""]

# List the unique ASVs at each locaton

S1_ASV <- unique(amp_S1$tax$OTU)
S2_ASV <- unique(amp_S2$tax$OTU)
N_ASV <- unique(amp_N$tax$OTU)


# Using eulerr, create euler diagrams for taxa and ASVs

add_quantity_line_break <- function(v) {
  tags <- v$children$canvas.grob$children$diagram.grob.1$children$tags$children

  tags <- do.call(grid::gList, lapply(tags, function(x) {
    x$children[[2]]$label <- sub(" \\%)", "%)", x$children[[2]]$label)
    before <- gsub(" .*$", "", x$children[[2]]$label)
    after <- gsub(".*\\(", "(", x$children[[2]]$label)
    full <- bquote(atop(NA, atop(atop(textstyle(.(before)), textstyle(italic(.(after)))), NA)))

    x$children[[2]]$label <- full
    x$children[[2]]$just <- NULL
    x
  }))

  v$children$canvas.grob$children$diagram.grob.1$children$tags <- tags

  return(v)
}

E_taxa <- euler(list(S1 = S1_taxa, S2 = S2_taxa, N = N_taxa), shape = "ellipse")
E_taxa_plot <- plot(E_taxa, quantities = list(type = c("counts", "percent"), cex = 1), legend = list(side = "right", nrow = 1, ncol = 3, cex = 1))
E_taxa_plot <- add_quantity_line_break(E_taxa_plot)

unique_taxa_location <- (E_taxa$original.values["S1"] +
  E_taxa$original.values["S2"] +
  E_taxa$original.values["N"]) / sum(E_taxa$original.values)


E_ASV <- euler(list(S1 = S1_ASV, S2 = S2_ASV, N = N_ASV), shape = "ellipse")
E_ASV_plot <- plot(E_ASV, quantities = list(type = c("counts", "percent"), cex = 1), legend = list(side = "right", nrow = 1, ncol = 3, cex = 1))
E_ASV_plot <- add_quantity_line_break(E_ASV_plot)

unique_ASV_location <- (E_ASV$original.values["S1"] +
  E_ASV$original.values["S2"] +
  E_ASV$original.values["N"]) / sum(E_ASV$original.values)

Across both taxa (n=r unique["Species"]``unique["Species"]) and ASVs (n=27038 processed_ASVs), all three locations had unique elements, suggesting that there was not complete mixing of eDNA across micro-habitats. 33.3% sprintf("%0.1f%%", unique_taxa_location * 100) of taxa and 79.1% sprintf("%0.1f%%", unique_ASV_location * 100) of ASVs were unique to one of the three locations (Figure 5). At both the taxa and ASV level, S1 had more unique elements than S2 and nearshore, although this effect was more pronounced across ASVs. ASVs and taxa shared across pairs of two locations were the least common elements.

euler <- (wrap_elements(panel = grid::textGrob("A)", hjust = 6, vjust = 1, gp = gpar(fontsize = 14, fontface = "bold"))) +
  wrap_elements(panel = grid::textGrob("B)", hjust = 6, vjust = 1, gp = gpar(fontsize = 14, fontface = "bold")))) /
  (wrap_elements(panel = E_taxa_plot$children$canvas.grob) + wrap_elements(panel = E_ASV_plot$children$canvas.grob)) /
  wrap_elements(panel = E_taxa_plot$children$legend.grob) + plot_layout(heights = c(1, 10, 1))

path <- "Figures/Figure5."
save_pdf_png(plot = euler, path = path, w = 169, h = 110, u = "mm", bg = "transparent", dpi = 600)

Figure 5: Area-proportional Euler diagrams showing the overlap in taxa (A) and ASVs (B) between locations S1, S2, and N.

# Function for creating indicator species heatmap from physeq data

indicator_heatmap <- function(physeq, asv) {
  # If analyzing by taxa, remove all unclassified ASVs and condense to species assignments
  if (asv == FALSE) {
    physeq <- tax_glom(physeq, taxrank = "Species", NArm = TRUE)
    tax_table(physeq) <- remove_NCBInums(tax_table(physeq))
    # tax_table(physeq) = gsub("\\_.*","",tax_table(physeq))
  } else {
    tax_table(physeq) <- remove_NCBInums(tax_table(physeq))
  }

  # Convert to presence/absence data
  otu_table(physeq)[otu_table(physeq) > 0] <- 1

  # Use indicspecies package to get top taxa
  indval <- multipatt(t(otu_table(physeq)), t(sample_data(physeq)[, "Location"]), control = how(nperm = 999), duleg = TRUE)
  top_N <- rownames(indval$sign %>% filter(p.value < 0.05) %>% filter(s.N == 1) %>% filter(stat > .7))
  top_S1 <- rownames(indval$sign %>% filter(p.value < 0.05) %>% filter(s.S1 == 1) %>% filter(stat > .7))
  top_S2 <- rownames(indval$sign %>% filter(p.value < 0.05) %>% filter(s.S2 == 1) %>% filter(stat > .7))

  indicators_all <- c(top_N, top_S1, top_S2)

  # Merge samples by site/time before plotting
  physeq_merge <- merge_samples(physeq, "SiteByTime")

  # Update metadata accordingly
  updated_metadata <- as.matrix(sample_data(physeq)[, c("SiteByTime", "Time", "Location", "T", "S")])
  updated_metadata <- as.data.frame(updated_metadata)
  updated_metadata$SiteByTime <- as.factor(updated_metadata$SiteByTime)
  updated_metadata$Time <- as.factor(updated_metadata$Time)
  updated_metadata$Location <- as.factor(updated_metadata$Location)
  updated_metadata <- distinct(updated_metadata)
  rownames(updated_metadata) <- updated_metadata[, "SiteByTime"]
  sample_data(physeq_merge) <- sample_data(updated_metadata)

  # Divide all merged samples by the number of samples they represent
  otu_table(physeq_merge) <- otu_table(physeq_merge) / 3

  otu_table(physeq_merge)[grepl("11:30", rownames(otu_table(physeq_merge)))] <- otu_table(physeq_merge)[grepl("11:30", rownames(otu_table(physeq_merge)))] / 3

  otu_table(physeq_merge)[grepl("14:00", rownames(otu_table(physeq_merge)))] <- otu_table(physeq_merge)[grepl("14:00", rownames(otu_table(physeq_merge)))] / 3

  otu_table(physeq_merge)[grepl("17:00", rownames(otu_table(physeq_merge)))] <- otu_table(physeq_merge)[grepl("17:00", rownames(otu_table(physeq_merge)))] / 3

  # Add labels to the indicator species for plotting
  indicators_only <- prune_taxa(colnames(otu_table(physeq_merge)) %in% indicators_all, physeq_merge)

  add_indicator_label <- function(X) {
    if (X %in% top_N) {
      return("N")
    } else if (X %in% top_S1) {
      return("S1")
    } else if (X %in% top_S2) {
      return("S2")
    } else {
      return("NA")
    }
  }

  indicators <- sapply(rownames(tax_table(indicators_only)), add_indicator_label)

  tax_table(indicators_only) <- cbind(tax_table(indicators_only), indicators)

  # Plot
  plot <- plot_heatmap(indicators_only,
    method = NULL,
    distance = NULL,
    trans = NULL,
    taxa.order = rev(indicators_all),
    taxa.label = "Species",
    sample.label = "Time",
    sample.order = "Time"
  ) +
    facet_grid(factor(indicators, levels = c("S1", "S2", "N")) ~ factor(Location, levels = c("S1", "S2", "N")),
      scales = "free",
      space = "free",
      switch = "y"
    ) +
    scale_fill_gradient(low = "black", high = "#B5D7E4") +
    theme_grey(base_size = 10) +
    theme(plot.title = element_text(hjust = 0.5)) +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
    theme(strip.placement = "outside") +
    theme(strip.text = element_text(size = rel(1))) +
    theme(legend.position = "bottom") +
    theme(panel.background = element_blank(), panel.grid.major = element_blank()) +
    labs(fill = "Proportion of Positive Detects", title = "Locations") +
    theme(
      axis.ticks = element_line(color = "black"),
      axis.text.y = element_text(color = "black", face = "italic"),
      axis.text.x = element_text(color = "black")
    ) +
    guides(fill = guide_colorbar(barwidth = 10))

  plot$scales$scales[[2]]$name <- "Indicators"

  return(list(
    plot = plot,
    top_N = top_N,
    top_S1 = top_S1,
    top_S2 = top_S2,
    all_indicators = indicators_all,
    full_output = indval
  ))
}

heatmap_taxa <- indicator_heatmap(data_physeq, FALSE)
heatmap_ASV <- indicator_heatmap(data_physeq, TRUE)

indicator_ASVs <- data_amp$tax[heatmap_ASV$all_indicators, ]
unassigned_indicator_ASVs <- indicator_ASVs %>%
  group_by(Kingdom) %>%
  summarise(n = n()) %>%
  mutate(freq = n / sum(n))


identify_new_indicators <- function(data_amp, top_ASVs, top_taxa) {
  top_ASVs <- data_amp$tax[top_ASVs, ]
  top_taxa <- data_amp$tax[top_taxa, ]

  top_ASVs <- top_ASVs$Species[top_ASVs$Species != ""]
  top_taxa <- top_taxa$Species

  top_ASVs <- remove_NCBInums(top_ASVs)
  top_taxa <- remove_NCBInums(top_taxa)
  return(setdiff(top_ASVs, top_taxa))
}

N_ASV_extra <- sort(identify_new_indicators(data_amp, heatmap_ASV$top_N, heatmap_taxa$top_N))
S1_ASV_extra <- sort(identify_new_indicators(data_amp, heatmap_ASV$top_S1, heatmap_taxa$top_S1))
S2_ASV_extra <- sort(identify_new_indicators(data_amp, heatmap_ASV$top_S2, heatmap_taxa$top_S2))

path <- "Figures/Figure6."
save_pdf_png(plot = heatmap_taxa$plot, path = path, w = 169, h = 100, u = "mm", bg = "transparent", dpi = 600)

Figure 6: Heat map of identified indicator taxa for all three locations over time, scaled by the number of replicates that detected that species.

Using an indicator species framework produced similar trends to the analysis of unique taxa and ASVs; all three locations had indicator taxa and ASVs, and S1 had more indicator elements than other locations. By taxa, we identified 11 length(heatmap_taxa$all_indicators) indicator taxa: 7 length(heatmap_taxa$top_S1) for S1, 3 length(heatmap_taxa$top_S2) for S2, and 1 length(heatmap_taxa$top_N) for nearshore (Figure 6). By ASVs, we found more indicators overall. We identified 99 length(heatmap_ASV$all_indicators) indicator ASVs: 72 length(heatmap_ASV$top_S1) for S1, 10 length(heatmap_ASV$top_S2) for S2, and 17 length(heatmap_ASV$top_N) for nearshore (Supplemental Figure X). Most (78.8% sprintf("%0.1f%%", unassigned_indicator_ASVs[unassigned_indicator_ASVs$Kingdom == "",]$freq * 100)) indicator ASVs had no taxonomic assignment, but the ones that did included several species beyond the original indicator taxa, at S1 (Bossiella frondescens, Bossiella plumosa, Corallina confusa, Leathesia difformis, Myrionema balticum, Scorpaenichthys marmoratus toString(S1_ASV_extra)), S2 (Bossiella frondescens, Smithora naiadum toString(S2_ASV_extra)), and nearshore (Halichondria panicea, Mazzaella splendens, Smithora naiadum, Strongylocentrotus purpuratus toString(N_ASV_extra)).

3.3. Community-Level Differences between Locations

# Calculate Jaccard dissimilarity matrices
jaccard_full <- vegdist(t(data_amp$abund), binary = TRUE, method = "jaccard")
jaccard_taxa <- vegdist(t(data_amp_taxa$abund), binary = TRUE, method = "jaccard")

set.seed(0)

# PERMANOVA on full data set
permutest(betadisper(jaccard_full, data_amp$metadata$Location))

## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##            Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)  
## Groups      2 0.002340 0.00117003 2.9267    999  0.061 .
## Residuals 156 0.062365 0.00039977                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(jaccard_full ~ Location + Time + Dereplicated_Sample_Name, data = data_amp$metadata)

# PERMANOVA on taxa data set
permutest(betadisper(jaccard_taxa, data_amp_taxa$metadata$Location))

## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##            Df  Sum Sq    Mean Sq      F N.Perm Pr(>F)
## Groups      2 0.00061 0.00030272 0.1364    999  0.894
## Residuals 156 0.34626 0.00221963

adonis2(jaccard_taxa ~ Location + Time + Dereplicated_Sample_Name, data = data_amp_taxa$metadata)

# Check the stress vs. dimension plot

# Function that performs a NMDS for 1-10 dimensions and plots the nr of dimensions vs the stress
NMDS.scree <- function(x) { # where x is the name of the data frame variable
  plot(rep(1, 10), replicate(10, metaMDS(x, autotransform = F, k = 1)$stress), xlim = c(1, 10), ylim = c(0, 0.30), xlab = "# of Dimensions", ylab = "Stress", main = "NMDS stress plot")
  for (i in 1:10) {
    points(rep(i + 1, 10), replicate(10, metaMDS(x, autotransform = F, k = i + 1)$stress))
  }
}

# Use NMDS.scree function to choose the optimal nr of dimensions
NMDS.scree(jaccard_full)

NMDS.scree(jaccard_taxa)

# Establish functions for PAM clustering algorithm

PAM_validation <- function(dist) {
  k_list <- c()
  avg_sil_width <- c(0)
  for (x in 1:(nrow(as.matrix(dist)) - 1)) {
    PAM <- pam(dist, k = x, diss = TRUE)
    avg_sil_width <- c(avg_sil_width, PAM$silinfo$avg.width)
    k_list <- c(k_list, x)
  }

  PAM_validation <- data.frame(k = k_list, sil_width = avg_sil_width, x_max = which.max(avg_sil_width))
  return(PAM_validation)
}

PAM_validation_viz <- function(PAM_validation) {
  plot <- ggplot(
    data = PAM_validation,
    aes(x = k, y = sil_width)
  ) +
    geom_line() +
    geom_point() +
    xlab("K") +
    ylab("Average Silhouette Width") +
    ggtitle("Cluster Validation using Average Silhouette Width") +
    geom_segment(
      aes(
        x = x_max[1],
        y = 0,
        xend = x_max[1],
        yend = sil_width[x_max[1]],
        color = "red"
      ),
      linetype = 2,
      show.legend = FALSE
    )
  return(plot)
}

PAM_silhouette_viz_3 <- function(PAM, sample_data) {
  set.seed(0)
  PAM_df <- data.frame(PAM$silinfo$widths)
  PAM_df <- cbind(PAM_df, X = row.names(PAM_df))
  PAM_df <- merge(PAM_df, sample_data, by = "X")
  PAM_df$cluster <- as.factor(PAM_df$cluster)
  PAM_df$random <- sample(1000, size = nrow(PAM_df), replace = FALSE)


  ggplot(
    data = PAM_df,
    aes(x = sil_width, y = reorder(X, random), fill = Location)
  ) +
    geom_point(aes(shape = Location, color = Time), size = 2, position = position_jitter(height = 20, width = 0)) +
    scale_shape_manual(values = c(16, 15, 17)) +
    facet_grid(cluster ~ ., scales = "free_y", space = "free_y") +
    viridis::scale_color_viridis(discrete = TRUE, begin = 0, end = .97, direction = -1) +
    xlab("Silhouette Width") +
    ylab("") +
    geom_vline(xintercept = PAM$silinfo$avg.width, linetype = "dashed", color = "red") +
    theme_classic(base_size = 10) +
    theme(axis.text.x = element_text(angle = 0)) +
    ggtitle("PAM Clustering Visualization") +
    theme(
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank(),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.background = element_blank(),
      panel.border = element_rect(colour = "black", size = 1, fill = NA)
    ) +
    scale_y_discrete(expand = expansion(add = 40))
}

NMDS_PAM <- function(amp, jaccard, dim, PAM_logical) {
  set.seed(0)
  NMDS <- metaMDS(jaccard, dim, trymax = 100, tidy = TRUE)

  NMDS_stress <- stressplot(NMDS)

  data.scores <- as.data.frame(scores(NMDS))
  data.scores$X <- rownames(data.scores)
  data.scores <- merge(data.scores, amp$metadata)

  annotations1 <- data.frame(
    xpos = c(Inf), ypos = c(-Inf),
    annotateText = c(paste("Stress:", round(NMDS$stress, 3))),
    hjustvar = c(1),
    vjustvar = c(-.5)
  )

  plot1 <- ggplot(data = data.scores, aes(x = NMDS1, y = NMDS2)) +
    theme_light() +
    geom_point(aes(shape = Location, color = Time), size = 2) +
    scale_shape_manual(values = c(16, 15, 17)) +
    theme(text = element_text(size = 10)) +
    ggtitle("A) nNMDS Ordination (Jaccard Dissimilarity)") +
    viridis::scale_color_viridis(discrete = TRUE, begin = 0, end = .97, direction = -1) +
    stat_ellipse(aes(group = Location), type = "t", linetype = 2, alpha = 1)


  plot1 <- plot1 +
    geom_text(data = annotations1, aes(x = xpos, y = ypos, hjust = hjustvar, vjust = vjustvar, label = annotateText)) +
    theme_classic(base_size = 10)

  if (PAM_logical == TRUE) {
    validated_k <- PAM_validation(jaccard)
    PAM_validation <- PAM_validation_viz(validated_k)

    PAM <- pam(jaccard, k = validated_k$x_max[1], diss = TRUE)
    plot2 <- PAM_silhouette_viz_3(PAM, amp$metadata)

    return(list(NMDS = plot1, PAM = plot2, NMDS_stress = NMDS_stress, PAM_validation = PAM_validation))
  }

  return(plot1)
}


NMDS_PAM_full <- NMDS_PAM(data_amp, jaccard_full, 2, TRUE)

NMDS_PAM_full$PAM_validation

NMDS_PAM_taxa <- NMDS_PAM(data_amp_taxa, jaccard_taxa, 2, TRUE)

NMDS_PAM_taxa$PAM_validation

NMDS_PAM_combined <- (NMDS_PAM_full$NMDS + theme(legend.position = "none") + ggtitle("A)")) +
  (NMDS_PAM_full$PAM + theme(legend.position = "none") + ggtitle("B)")) +
  (NMDS_PAM_taxa$NMDS + theme(legend.position = "none") + ggtitle("C)")) +
  (NMDS_PAM_taxa$PAM + ggtitle("D)")) +
  plot_layout(guides = "collect", widths = c(1, 1))

path <- "Figures/Figure7."
save_pdf_png(plot = NMDS_PAM_combined, path = path, w = 169, h = 169, u = "mm", bg = "transparent", dpi = 600)

Community composition also differed across locations, as supported by multiple analyses. Variation in community composition was significantly explained by location, even when accounting for variation due to time of sampling and variation between replicates (ASVs: PERMANOVA p < 0.001, betadisper p > 0.05; Taxa: PERMANOVA p < 0.001, betadisper p > 0.05). As visualized using NMDS ordination, S1 and S2 samples collected around low tide were most dissimilar from the offshore samples (Figures 7A, 7C). These patterns persisted in the optimal clusters identified using the partitioning among medoids (PAM) algorithm, without assuming a priori that location drives the clusters (Figure 7B, 7D). When analyzed by ASV, PAM identified three optimal clusters, primarily composed of: 1) nearshore samples with additional samples from S1 and S2 at early time points, 2) the remaining S1 samples, and 3) the remaining S2 samples (Figure 7B). When analyzing by taxa, PAM additionally differentiated between samples from the three locations collected at early time points. Thus, multiple analyses with different initial assumptions all demonstrated differentiable eDNA signals across micro-habitats.

Figure 7: : Ordination analysis (NMDS using Jaccard Distance; A, C) and cluster analysis (partitioning among medoids with number of clusters validated to maximize average silhouette width; B,D). Panels A & B are analyzed by ASV, and panels C & D are analyzed by taxa. 95% confidence ellipses for the centroids of each location based on a multivariate t-distribution are shown in the NMDS ordinations with black dashed lines, and average silhouette width is shown in cluster diagrams with a red dashed line.

distmat <- as.matrix(jaccard_full)

# Calculate Pairwise Dissimilarity Values Within and Across Sites

Jaccard <- c()
Combo <- c()
Time <- c()
summary <- data.frame(Jaccard, Combo, Time)

one <- c("S1", "S2", "N", "S1", "S2", "S1")
one_names <- c("S1", "S2", "N", "S1", "S2", "S1")
two <- c("S1", "S2", "N", "N", "N", "S2")
two_names <- c("S1", "S2", "N", "N", "N", "S2")

for (x in levels(data_amp$metadata$Time)) {
  subset <- filter(data_amp$metadata, data_amp$metadata$Time == x)

  for (i in 1:length(one)) {
    name <- paste(one_names[i], two_names[i], sep = ":")
    data <- distmat[subset$X[subset$Location == one[i]], subset$X[subset$Location == two[i]]]
    Jaccard <- vector()

    if (!(0 %in% dim(data))) {
      for (row in 1:nrow(data)) {
        for (col in 1:ncol(data)) {
          if (col > row) {
            Jaccard <- c(Jaccard, data[row, col])
          }
        }
      }
    }
    Combo <- rep(name, length(Jaccard))
    Time <- rep(x, length(Jaccard))

    summary <- rbind(summary, data.frame(Jaccard, Combo, Time))
  }
}


## PLOTTING

# Define palette for colors and fills
colors <- turbo(8)

palette <- c(
  "Within Site" = "#30123BFF",
  "S1:N" = "#1AE4B6FF",
  "S2:N" = "#FABA39FF",
  "S1:S2" = "#7A0403FF"
)

fill_palette <- c(
  "Within Site" = "#30123B20",
  "S1:N" = "#1AE4B620",
  "S2:N" = "#FABA3920",
  "S1:S2" = "#7A040320"
)


# WITHIN SITE COMBINED

summary$Combo2 <- summary$Combo
summary$Combo2[summary$Combo2 == "S1:S1"] <- "Within Site"
summary$Combo2[summary$Combo2 == "S2:S2"] <- "Within Site"
summary$Combo2[summary$Combo2 == "N:N"] <- "Within Site"

stats <- summary %>%
  filter(Combo2 == "Within Site") %>%
  summarize(mean = mean(Jaccard), max = max(Jaccard))


g2.5 <- ggplot(summary, aes(
  x = Time,
  y = Jaccard,
  color = factor(Combo2, level = c("Within Site", "S1:N", "S2:N", "S1:S2")),
  fill = factor(Combo2, level = c("Within Site", "S1:N", "S2:N", "S1:S2"))
)) +
  geom_boxplot(position = position_dodge(.5, preserve = "single")) +
  geom_point(position = position_jitterdodge(dodge.width = .5), alpha = 0.2) +
  scale_color_manual(values = palette, name = "Site Pairings") +
  scale_fill_manual(values = fill_palette, name = "Site Pairings") +
  geom_hline(yintercept = stats$mean, color = palette["Within Site"]) +
  ylab("Jaccard Dissimilarity") +
  coord_cartesian(expand = 0, clip = "off") +
  geom_text(
    aes(max(Time),
      stats$mean,
      label = "Mean Dissimilarity\n  Within Sites",
      hjust = -.2
    ),
    color = palette["Within Site"],
    show.legend = FALSE,
    check_overlap = TRUE
  ) +
  theme_classic(base_size = 12) +
  theme(
    axis.ticks = element_line(color = "black"),
    axis.text.y = element_text(color = "black"),
    axis.text.x = element_text(color = "black")
  ) +
  theme(legend.position = "right", legend.justification = "top")

path <- "Figures/Figure8."
save_pdf_png(plot = g2.5, path = path, w = 169, h = 120, u = "mm", bg = "transparent", dpi = 600)

The extent of the Jaccard dissimilarity between samples from different locations varied over time. As shown in Figure 8, except at the first time point sampled, the dissimilarity between samples across sites was always greater than the mean dissimilarity between samples within the same site, and increasingly so across the period sampled. At two time points—15:30 and 16:00—after low tide but before water was moving between all locations again, all dissimilarity values across sites were greater than the maximum dissimilarity recorded between two samples from the same site.

Figure 8: Box plots over time showing the pairwise Jaccard dissimilarity values between all unique pairs of replicates from the same sites (Within Site) and all unique pairs of samples across two sites (S1:O, S:O, S1:S2). The solid line depicts the mean Jaccard dissimilarity value from within site pairings.

3.4. Ecological Significance

### Ecological Significance###

# NOT DONE YET

S1_unique <- setdiff(S1_taxa, union(S2_taxa, N_taxa))
S2_unique <- setdiff(S2_taxa, union(S1_taxa, N_taxa))
N_unique <- setdiff(N_taxa, union(S2_taxa, S1_taxa))


S1 <- data.frame(list(Species = remove_NCBInums(S1_unique), Unique = "S1"))
S2 <- data.frame(list(Species = remove_NCBInums(S2_unique), Unique = "S2"))
N <- data.frame(list(Species = remove_NCBInums(N_unique), Unique = "N"))
unique <- rbind(S1, S2, N)

S1_i <- data.frame(list(Species = remove_NCBInums(data_amp$tax$Species[rownames(data_amp$tax) %in% heatmap_taxa$top_S1]), Indicator = "S1"))
S2_i <- data.frame(list(Species = remove_NCBInums(data_amp$tax$Species[rownames(data_amp$tax) %in% heatmap_taxa$top_S2]), Indicator = "S2"))
N_i <- data.frame(list(Species = remove_NCBInums(data_amp$tax$Species[rownames(data_amp$tax) %in% heatmap_taxa$top_N]), Indicator = "N"))
indicators <- rbind(S1_i, S2_i, N_i)

df_list <- list(unique, indicators)
full_taxa_list <- df_list %>% reduce(full_join, by = "Species")





synonyms <- c()
for (s in full_taxa_list$Species) {
  print(s)
  worms_ID <- get_wormsid(s)
  if (is.na(worms_ID)) {
    syns <- NA
  } else {
    syns <- synonyms(worms_ID, db = "worms")
    syns <- toString(syns[[1]]$scientificname)
  }
  synonyms <- rbind(synonyms, syns)
}


full_taxa_list <- cbind(full_taxa_list, synonyms)

# Remove anything in parentheticals
full_taxa_list$synonyms <- gsub("\\s*\\([^\\)]+\\)", "", full_taxa_list$synonyms)

# Remove anything beyond two words in each comma separated section
for (syn in 1:length(full_taxa_list$synonyms)) {
  print(syn)
  all <- unlist(strsplit(full_taxa_list$synonyms[syn], ","))
  if (!is.na(all[1])) {
    print(all)
    for (i in 1:length(all)) {
      print(all[i])
      all[i] <- str_trim(all[i])
      all[i] <- sub("^(\\S*\\s+\\S+).*", "\\1", all[i])
      print(all[i])
    }
  }
  syns <- toString(unique(all))
  full_taxa_list$synonyms[syn] <- syns
}


write.csv(full_taxa_list, "Analysis Products/BPT_Analysis.csv", row.names = FALSE)


# Read in needed text of Between Pacific Tides (file not given on GitHub due to copyright)
txt <- pdf_text("Data/BetweenPacificTides_TextForSearching.pdf")
toc <- pdf_text("Data/BetweenPacificTides_TOC.pdf")

# first remove hyphenated situations "-\n"
# then remove all line breaks "\n"
txt <- gsub("-\n", "", txt)
txt <- gsub("\n", " ", txt)

# Clean up TOC
toc <- strsplit(toc, "\n")
toc <- c(toc[[1]], toc[[2]], toc[[3]])
page_nums <- grepl("[1234567890,]$", toc)
toc <- toc[page_nums]

# Fix the couple of strings that are split over multiple lines
toc_cleaned <- c(toc[1:14], paste0(toc[15], toc[16]), toc[17:21], paste0(toc[22], toc[23]), toc[24:53])
toc_cleaned <- as.data.frame(toc_cleaned)
toc_cleaned <- extract(toc_cleaned, toc_cleaned, into = c("text", "page"), "(.*)\\s+([^ ]+)$")
toc_cleaned$page <- as.numeric(toc_cleaned$page)


toc_cleaned_max <- vector()
for (x in 1:(nrow(toc_cleaned) - 1)) {
  print(x)
  if (toc_cleaned$page[x] == toc_cleaned$page[x + 1]) {
    toc_cleaned_max[x] <- toc_cleaned$page[x]
  } else {
    toc_cleaned_max[x] <- toc_cleaned$page[x + 1] - 1
  }
}
toc_cleaned_max[nrow(toc_cleaned)] <- NA

toc_cleaned$max_page <- toc_cleaned_max

# remove everything not referring to rocky shores
toc_cleaned <- toc_cleaned[1:17, ]

# Remove just rows with zone
toc_zones <- toc_cleaned %>% filter(grepl("Zone", text))
toc_zones$text <- gsub("\\..*", "", toc_zones$text)
toc_zones$text <- trimws(toc_zones$text)


output <- c()
for (x in 1:length(full_taxa_list$Species)) {
  species <- full_taxa_list$Species[x]
  synonyms <- full_taxa_list$synonyms[x]
  synonyms <- unlist(strsplit(synonyms, ", "))
  all <- c(species, synonyms)
  all_abbreviated <- sub("[a-z]* ", ". ", all)
  all_to_search <- c(all, all_abbreviated)

  all_pages <- vector()
  for (y in all_to_search) {
    print(y)
    pages <- grep(y, txt, ignore.case = FALSE, fixed = TRUE)
    all_pages <- c(all_pages, pages)
    print(pages)
  }
  all_pages <- unique(all_pages)
  output[x] <- list(all_pages)
}

full_taxa_list$pages <- output
full_taxa_list$pages[sapply(full_taxa_list$pages, function(x) length(x) == 0)] <- NA



zones <- unique(toc_zones$text)
df <- data.frame(matrix(ncol = length(zones), nrow = 0))
colnames(df) <- zones

name <- c()
zone <- c()
page <- c()

for (x in 1:length(full_taxa_list$pages)) {
  list <- full_taxa_list$pages[[x]]
  if (!is.na(list[1])) {
    for (y in 1:length(list)) {
      num <- list[y]
      for (z in 1:length(toc_zones$text)) {
        min <- toc_zones$page[z]
        max <- toc_zones$max[z]
        if (num <= max && num >= min) {
          name <- append(name, full_taxa_list$Species[x])
          zone <- append(zone, toc_zones$text[z])
          page <- append(page, num)
        } else {
        }
      }
    }
  }
}

final_1 <- as.data.frame(cbind(name, zone, page))

# Split "Zones 1 and 2" into two entries
to_split <- final_1 %>%
  filter(zone == "Zones 1 and 2") %>%
  mutate(zone = "Zone 1")
final_2 <- final_1 %>% mutate(zone = replace(zone, zone == "Zones 1 and 2", "Zone 2"))
final_2 <- rbind(final_2, to_split)


# Split "Zones 2 and 3" into two entries
to_split <- final_1 %>%
  filter(zone == "Zones 2 and 3") %>%
  mutate(zone = "Zone 2")
final_2 <- final_2 %>% mutate(zone = replace(zone, zone == "Zones 2 and 3", "Zone 3"))
final_2 <- rbind(final_2, to_split)


condense <- final_2 %>% pivot_wider(names_from = zone, values_from = page, values_fn = unique(list()))


condense <- left_join(condense, full_taxa_list[, c("Species", "Unique", "Indicator")], by = join_by(name == Species))

condense <- condense %>% mutate(Site = coalesce(Unique, Indicator))
condense <- condense %>% relocate(Site, Unique, Indicator, name, "Zone 4", "Zone 3", "Zone 2", "Zone 1")

condense <- condense %>% rename(Zone_1 = "Zone 1", Zone_2 = "Zone 2", Zone_3 = "Zone 3", Zone_4 = "Zone 4")



condense_no_pages <- condense
condense_no_pages$Zone_1 <- !sapply(condense_no_pages$Zone_1, is.null)
condense_no_pages$Zone_2 <- !sapply(condense_no_pages$Zone_2, is.null)
condense_no_pages$Zone_3 <- !sapply(condense_no_pages$Zone_3, is.null)
condense_no_pages$Zone_4 <- !sapply(condense_no_pages$Zone_4, is.null)


# Adjust Based on Manual Checking

# remove "Proboscidactyla flavicirrata" because search identified the wrong abbreviated synonym
condense_no_pages <- condense_no_pages %>% filter(name != "Proboscidactyla flavicirrata")

# add "Laminaria setchellii" because it was found manually and only appears in a diagram (so not in text search)
condense_no_pages <- condense_no_pages %>% add_row(Site = "N", Unique = "N", Indicator = NA, name = "Laminaria setchellii", Zone_4 = TRUE, Zone_3 = FALSE, Zone_2 = FALSE, Zone_1 = FALSE)

# add "Egregia menziesii" because it was found manually through the index, but only appears by genus in the text
condense_no_pages <- condense_no_pages %>% add_row(Site = "S2", Unique = "S2", Indicator = NA, name = "Egregia menziesii", Zone_4 = TRUE, Zone_3 = TRUE, Zone_2 = TRUE, Zone_1 = FALSE)

# add "Pandalus danae" because it was found manually through the index, but only appears by genus in the text
condense_no_pages <- condense_no_pages %>% add_row(Site = "S2", Unique = "S2", Indicator = NA, name = "Pandalus danae", Zone_4 = TRUE, Zone_3 = FALSE, Zone_2 = FALSE, Zone_1 = FALSE)

# Adjust "Ectopleura marina" to include Zone_3 (found manually because mentioned on that page with genus only, but linked in the index)
condense_no_pages <- condense_no_pages %>% mutate(Zone_3 = replace(Zone_3, name == "Ectopleura marina", TRUE))

# Adjust "Halosydna brevisetosa" to include Zone_4 (found manually because mentioned on that page with genus only, but linked in the index)
condense_no_pages <- condense_no_pages %>% mutate(Zone_4 = replace(Zone_4, name == "Halosydna brevisetosa", TRUE))



condense_no_pages <- condense_no_pages %>% arrange(name)
condense_no_pages <- condense_no_pages %>% arrange(desc(Zone_1), desc(Zone_2), desc(Zone_3), desc(Zone_4))
condense_no_pages <- condense_no_pages %>% arrange(factor(Site, levels = c("S1", "S2", "N")))







contingency_table <- condense_no_pages %>%
  group_by(Site) %>%
  summarise(High = sum(Zone_1 == TRUE | Zone_2 == TRUE | Zone_3 == TRUE), Low = sum(Zone_1 == FALSE & Zone_2 == FALSE & Zone_3 == FALSE))
contingency_table <- contingency_table %>%
  select(-Site) %>%
  as.matrix()
rownames(contingency_table) <- c("N", "S1", "S2")


set.seed(1)

test <- chisq.test(contingency_table, correct = FALSE, simulate.p.value = TRUE)

test_post_hoc <- chisq.posthoc.test(contingency_table, method = "bonferroni", simulate.p.value = TRUE)


for_plotting <- condense_no_pages
for_plotting <- for_plotting %>% mutate(Unique = replace(Unique, !is.na(Unique), "X"), Unique = replace(Unique, is.na(Unique), ""), Indicator = replace(Indicator, !is.na(Indicator), "X"), Indicator = replace(Indicator, is.na(Indicator), ""))




for_plotting <- for_plotting %>%
  gt(groupname_col = "Site", id = "mygt") %>%
  tab_spanner(
    label = "Intertidal Zones",
    columns = c(
      Zone_4, Zone_3, Zone_2, Zone_1
    )
  ) %>%
  tab_options(row_group.as_column = TRUE) %>%
  gt_color_rows("Zone_4", palette = c("white", "#D1E7C7")) %>%
  gt_color_rows("Zone_3", palette = c("white", "#BBBED1")) %>%
  gt_color_rows("Zone_2", palette = c("white", "#BBBED1")) %>%
  gt_color_rows("Zone_1", palette = c("white", "#BBBED1")) %>%
  tab_style(
    style = cell_text(style = "italic"),
    locations = cells_body(
      columns = name
    )
  ) %>%
  tab_style(
    style = cell_text(color = "#FFFFFF00", size = "x-small"),
    locations = cells_body(
      columns = c("Zone_4", "Zone_3", "Zone_2", "Zone_1")
    )
  ) %>%
  tab_style(
    style = cell_text(size = px(12)),
    locations = cells_body(
      columns = c("Unique", "Indicator", "name")
    )
  ) %>%
  cols_align(
    align = "center",
    columns = c("Unique", "Indicator")
  ) %>%
  cols_label(
    Zone_1 = "1",
    Zone_2 = "2",
    Zone_3 = "3",
    Zone_4 = "4",
    name = "Species",
    Unique = "Uniq.\nTo",
    Indicator = "Indic.\nOf"
  ) %>%
  tab_style(
    style = cell_text(v_align = "middle"),
    locations = cells_row_groups()
  ) %>%
  tab_style(
    style = cell_borders(
      sides = c("bottom"),
      color = "black",
      weight = px(3)
    ),
    locations = list(cells_column_spanners(), cells_column_labels())
  ) %>%
  tab_style(
    style = cell_borders(
      sides = "all",
      color = "black",
      weight = px(3)
    ),
    locations = list(cells_row_groups())
  ) %>%
  tab_style(
    style = cell_borders(
      sides = "all",
      color = "black"
    ),
    locations = cells_body()
  ) %>%
  tab_style(
    style = cell_borders(
      sides = "bottom",
      color = "black",
      weight = px(3)
    ),
    locations = cells_body(rows = c(18, 31, 41))
  ) %>%
  tab_style(
    style = cell_borders(
      sides = "right",
      color = "black",
      weight = px(3)
    ),
    locations = cells_body(columns = 8)
  ) %>%
  tab_options(data_row.padding = px(1)) %>%
  data_color(
    columns = Unique,
    rows = Unique == "X",
    palette = "lightgray",
  ) %>%
  data_color(
    columns = Indicator,
    rows = Indicator == "X",
    palette = "lightgray",
  ) %>%
  tab_options(
    table.border.top.style = "hidden"
  ) %>%
  tab_style(
    style = cell_text(align = "center"),
    locations = cells_column_labels(columns = name)
  )



path <- "Figures/Figure9.png"
gtsave(for_plotting, filename = path)

Identified unique and indicator taxa reflected known ecological differences between locations. Only a subset of unique and indicator taxa were described in Between Pacific Tides: 18.8% (9/48) of nearshore taxa, 32.1% (17/53) of S1 taxa, and 31.1% (14/45) of S2. However, across the described subset, more taxa from S1 were categorized to high and middle intertidal zones than S2 and N (Figure 9; not made in R). The proportion of taxa from the high and middle intertidal varied significantly across locations (χ2 =11.67, p<0.05), matching the environmental characteristics of the sites.

Figure 9: Chart habitat information about unique and indicator species present in Between Pacific Tides

Supporting Information

### Supporting Figure 1###

path <- "Figures/FigureS1."
save_pdf_png(plot = heatmap_ASV$plot, path = path, w = 169, h = 300, u = "mm", bg = "transparent", dpi = 600)

Figure S1: Heat map of identified indicator ASVs for all three locations over time, scaled by the number of replicates that detected that species.

## Supporting Table 1###

overlap_table

write.csv(overlap_table, "Analysis Products/TableS1.csv", row.names = FALSE)

By analyzing occurrence data from the Global Biodiversity Information Facility (on 2023-12-15 Sys.Date()), we confirmed that 100 percentages$count[6] of r sum(percentages$count)``sum(percentages$count) identified species (24.1% sprintf("%0.1f%%", percentages$per[6] * 100)) have occurrence records at Pillar Point, specifically. An additional 193 percentages$count[5] identified species (46.5% sprintf("%0.1f%%", percentages$per[5] * 100)) have occurrence records in the California Current System (CCS) and known ranges that encompass Pillar Point. Of the remaining species, some only fulfilled one of two criteria; 11 percentages$count[4] identified species (2.7% sprintf("%0.1f%%", percentages$per[4] * 100)) had occurrence records in the CCS but known ranges that did not encompass Pillar Point and 37 percentages$count[3] identified species (8.9% sprintf("%0.1f%%", percentages$per[3] * 100)) had ranges that encompass Pillar Point but no occurrence records in the CCS. Additionally, 39 percentages$count[1] identified species (9.4% sprintf("%0.1f%%", percentages$per[1] * 100)) had neither occurrence records in the CCS nor known ranges that encompassed Pillar Point. Finally, 15 percentages$count[2] identified species (3.6% sprintf("%0.1f%%", percentages$per[2] * 100)) lacked sufficient occurrence records in GBIF to make any designation, and 20 percentages$count[7] identified species (4.8% sprintf("%0.1f%%", percentages$per[7] * 100)) did not have species-level records in GBIF. Full details by species can be found in Table S2.

## Supporting Table 2###

GBIF_summary

write.csv(GBIF_summary, "Analysis Products/TableS2.csv", row.names = TRUE)

Code for Making Versions of Figures for PowerPoint

Figure 3

### Version of Figure 3 optimized for PowerPoint slides/animations###

graph <- sample_schema

if (!dir.exists("Figures/PowerPoint")) {
  dir.create("Figures/PowerPoint", recursive = TRUE)
}

filename <- paste0("sample_schema_full.pdf")
path <- paste0("Figures/PowerPoint/", filename)
pdf(path, width = 13, height = 4.3)
print(graph)
dev.off()

graph$layers[[1]] <- NULL

alpha <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

for (i in 1:12) {
  alpha[i] <- 1
  graph <- graph +
    scale_color_viridis_d(direction = -1, alpha = alpha) +
    scale_fill_viridis_d(direction = -1, alpha = alpha)
  filename <- paste0("sample_schema", i, ".pdf")
  path <- paste0("Figures/PowerPoint/", filename)
  pdf(path, width = 13, height = 4.3)
  print(graph)
  dev.off()
}

Figure 4

### Version of Figure 4 optimized for PowerPoint slides/animations###

p <- ggtree::ggtree(pt_physeq, layout = "rectangular")

ordered_names <- ggtree::get_taxa_name(p)
ordered_names <- unique(ordered_names)
ordered_names[match(rownames(tax_table(pt_physeq)), ordered_names)] <- tax_table(pt_physeq)[, "Phylum"]
ordered_names <- unique(ordered_names)

names(colors_ordered) <- ordered_names

p <- p + geom_fruit(geom = geom_tile, mapping = aes(fill = Phylum), width = 3, show.legend = FALSE) + scale_fill_manual(values = colors_ordered) + scale_color_manual(values = colors_ordered)

p <- p + ggtree::geom_hilight(data = phyla_nodes, mapping = aes(node = node, fill = node_phylum), show.legend = FALSE)

p <- p +
  theme(
    panel.background = element_rect(fill = "transparent"),
    plot.background = element_rect(fill = "transparent", color = NA),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.background = element_rect(fill = "transparent"),
    legend.box.background = element_rect(fill = "transparent")
  )

ggsave(
  "Figures/PowerPoint/phylo.pdf",
  plot = p,
  width = 4.49,
  height = 7.06,
  dpi = 300,
  bg = "transparent"
)

Figure 6

### Version of Figure 6 optimized for PowerPoint slides/animations###

graph <- heatmap_taxa$plot

path <- "Figures/PowerPoint/indicator_figure.png"
ggsave(filename = path, plot = graph, width = 13, height = 6.69)

Figure 7

### Version of Figure 7 optimized for PowerPoint slides/animations###

alpha <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
graph <- NMDS_PAM_full$NMDS +
  ggtitle("NMDS Ordination (Jaccard Dissimilarity)") +
  theme_classic(base_size = 18)

filename <- paste("NMDS_figure_full.pdf")
path <- paste("Figures/PowerPoint/", filename)
pdf(path, width = 11.97, height = 6.14)
print(graph)
dev.off()

graph$layers[[2]] <- NULL

graph <- graph +
  stat_ellipse(aes(group = Location), type = "t", linetype = 2, alpha = 0)

for (i in 1:12) {
  alpha[i] <- 1
  graph <- graph +
    viridis::scale_color_viridis(
      discrete = TRUE,
      begin = 0,
      end = .97,
      direction = -1,
      alpha = alpha
    )
  filename <- paste("NMDS_figure", i, ".pdf")
  path <- paste("Figures/PowerPoint/", filename)
  pdf(path, width = 11.97, height = 6.14)
  print(graph)
  dev.off()
  alpha[i] <- .2
}

Session Info

sessionInfo()

## R version 4.3.1 (2023-06-16)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.0.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices datasets  utils     methods  
## [8] base     
## 
## other attached packages:
##  [1] gtExtras_0.5.0           gt_0.10.0                pdftools_3.4.0          
##  [4] ampvis2_2.8.6            ggtreeExtra_1.12.0       ggtree_3.10.0           
##  [7] phyloseq_1.46.0          styler_1.10.2            devtools_2.4.5          
## [10] usethis_2.2.2            chunkhooks_0.0.1         ggh4x_0.2.6             
## [13] ggbeeswarm_0.7.2         mapview_2.11.2           rgbif_3.7.8             
## [16] mregions_0.1.8           spocc_1.2.2              pals_1.8                
## [19] Polychrome_1.5.1         taxize_0.9.100           chisq.posthoc.test_0.1.2
## [22] vegan_2.6-4              lattice_0.22-5           cluster_2.1.6           
## [25] indicspecies_1.7.14      permute_0.9-7            eulerr_7.0.0            
## [28] patchwork_1.1.3          sf_1.0-14                lubridate_1.9.3         
## [31] forcats_1.0.0            stringr_1.5.1            dplyr_1.1.4             
## [34] purrr_1.0.2              readr_2.1.4              tidyr_1.3.0             
## [37] tibble_3.2.1             ggplot2_3.4.4            tidyverse_2.0.0         
## 
## loaded via a namespace (and not attached):
##   [1] fs_1.6.3                bitops_1.0-7            bold_1.3.0             
##   [4] fontawesome_0.5.2       oai_0.4.0               httr_1.4.7             
##   [7] RColorBrewer_1.1-3      worrms_0.4.3            profvis_0.3.8          
##  [10] tools_4.3.1             utf8_1.2.4              R6_2.5.1               
##  [13] lazyeval_0.2.2          mgcv_1.9-0              measurements_1.5.1     
##  [16] rhdf5filters_1.14.1     urlchecker_1.0.1        withr_2.5.2            
##  [19] sp_2.1-2                gridExtra_2.3           prettyunits_1.2.0      
##  [22] leaflet_2.2.1           textshaping_0.3.7       leafem_0.2.3           
##  [25] cli_3.6.1               Biobase_2.62.0          labeling_0.4.3         
##  [28] sass_0.4.8              proxy_0.4-27            askpass_1.2.0          
##  [31] systemfonts_1.0.5       yulab.utils_0.1.0       R.utils_2.12.3         
##  [34] dichromat_2.0-0.1       sessioninfo_1.2.2       maps_3.4.1.1           
##  [37] rstudioapi_0.15.0       httpcode_0.3.0          generics_0.1.3         
##  [40] gridGraphics_0.5-1      crosstalk_1.2.1         Matrix_1.6-4           
##  [43] biomformat_1.30.0       fansi_1.0.5             S4Vectors_0.40.2       
##  [46] R.methodsS3_1.8.2       terra_1.7-55            lifecycle_1.0.4        
##  [49] whisker_0.4.1           scatterplot3d_0.3-44    yaml_2.3.7             
##  [52] rvertnet_0.8.2          rhdf5_2.46.1            paletteer_1.5.0        
##  [55] promises_1.2.1          crayon_1.5.2            miniUI_0.1.1.1         
##  [58] conditionz_0.1.0        chromote_0.1.2          geojson_0.3.5          
##  [61] mapproj_1.2.11          pillar_1.9.0            knitr_1.45             
##  [64] codetools_0.2-19        fastmatch_1.1-4         wk_0.9.1               
##  [67] glue_1.6.2              rebird_1.3.0            V8_4.4.1               
##  [70] ggfun_0.1.3             qpdf_1.3.2              data.table_1.14.8      
##  [73] remotes_2.4.2.1         urltools_1.7.3          vctrs_0.6.5            
##  [76] png_0.1-8               treeio_1.26.0           gtable_0.3.4           
##  [79] assertthat_0.2.1        rematch2_2.1.2          cachem_1.0.8           
##  [82] xfun_0.41               mime_0.12               survival_3.5-7         
##  [85] iterators_1.0.14        units_0.8-5             ellipsis_0.3.2         
##  [88] nlme_3.1-164            satellite_1.0.4         GenomeInfoDb_1.38.1    
##  [91] R.cache_0.16.0          bslib_0.6.1             vipor_0.4.5            
##  [94] KernSmooth_2.23-22      colorspace_2.1-0        BiocGenerics_0.48.1    
##  [97] DBI_1.1.3               raster_3.6-26           ade4_1.7-22            
## [100] phangorn_2.11.1         tidyselect_1.2.0        processx_3.8.2         
## [103] compiler_4.3.1          curl_5.1.0              geojsonsf_2.0.3        
## [106] xml2_1.3.6              plotly_4.10.3           triebeard_0.4.1        
## [109] scales_1.3.0            classInt_0.4-10         quadprog_1.5-8         
## [112] callr_3.7.3             digest_0.6.33           rmarkdown_2.25         
## [115] XVector_0.42.0          htmltools_0.5.7         pkgconfig_2.0.3        
## [118] base64enc_0.1-3         highr_0.10              fastmap_1.1.1          
## [121] rlang_1.1.2             htmlwidgets_1.6.4       shiny_1.8.0            
## [124] farver_2.1.1            jquerylib_0.1.4         zoo_1.8-12             
## [127] jsonlite_1.8.8          R.oo_1.25.0             RCurl_1.98-1.13        
## [130] magrittr_2.0.3          GenomeInfoDbData_1.2.11 ggplotify_0.1.2        
## [133] s2_1.1.4                Rhdf5lib_1.24.0         munsell_0.5.0          
## [136] Rcpp_1.0.11             viridis_0.6.4           ape_5.7-1              
## [139] ggnewscale_0.4.9        stringi_1.8.2           zlibbioc_1.48.0        
## [142] MASS_7.3-60             jqr_1.3.3               plyr_1.8.9             
## [145] pkgbuild_1.4.2          parallel_4.3.1          ggrepel_0.9.4          
## [148] geojsonio_0.11.3        Biostrings_2.70.1       splines_4.3.1          
## [151] multtest_2.58.0         hms_1.1.3               polylabelr_0.2.0       
## [154] ps_1.7.5                igraph_1.5.1            uuid_1.1-1             
## [157] reshape2_1.4.4          stats4_4.3.1            pkgload_1.3.3          
## [160] crul_1.4.0              ridigbio_0.3.7          evaluate_0.23          
## [163] leaflet.providers_2.0.0 renv_1.0.3              BiocManager_1.30.22    
## [166] tzdb_0.4.0              foreach_1.5.2           httpuv_1.6.13          
## [169] webshot2_0.1.1          polyclip_1.10-6         xtable_1.8-4           
## [172] e1071_1.7-14            tidytree_0.4.5          later_1.3.2            
## [175] viridisLite_0.4.2       class_7.3-22            ragg_1.2.6             
## [178] websocket_1.4.1         aplot_0.2.2             memoise_2.0.1          
## [181] beeswarm_0.4.0          IRanges_2.36.0          timechange_0.2.0